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Abstract 


The  project  New  Techniques  in  Time-Frequency  Analysis:  Adaptive  Band,  Ultra-Wide  Band  and  Multi-Rate 
Signal  Processing  led  to  the  development  of  new  techniques  and  theories  in  the  analysis  of  signals.  These 
techniques  and  theories  were  extensions  of  known  techniques  -  sampling,  Fourier,  Gabor  and  wavelet 
analysis,  and  new  approaches  to  analysis  -  using  combinations  of  analysis,  geometry,  group  theory,  and 
number  theory.  Every  target  item  in  the  original  statement  of  objectives  was  achieved,  and 
several  new  areas  of  research  were  open  up.  There  were  four  main  areas  of  study. 

Techniques  were  developed  to  deal  with  classes  of  signals  for  which  the  known  techniques  have  lim¬ 
itations,  e.g.,  adaptive  frequency  band  (AFB)  and  ultra- wide  band  (UWB)  signals.  These  signals  are 
of  interest  to  Air  Force  communications  systems,  but  the  theory  and  techniques  we  have  developed  and 
propose  to  continue  developing  fit  in  the  context  of  expanding  the  general  theory.  Our  techniques  in¬ 
volve  signal  segmentation,  projection,  expansion  in  bases,  transform  analysis,  and  multi-rate  methods. 
Adaptive  frequency  band  (AFB)  and  ultra- wide  band  (UWB)  systems  require  either  rapidly  changing  or 
very  high  sampling  rates.  Conventional  analog-to-digital  devices  have  non-adaptive  and  limited  dynamic 
range.  We  investigate  AFB  and  UWB  sampling  via  a  basis  projection  method.  The  method  decomposes 
the  signal  into  a  basis  over  time  segments  via  a  continuous-time  inner  product  operation  and  then  sam¬ 
ples  the  basis  coefficients  in  parallel.  The  signal  may  then  be  reconstructed  from  the  basis  coefficients  to 
recover  the  signal  in  the  time  domain.  We  develop  the  procedure  of  this  method,  analyze  various  meth¬ 
ods  for  signal  segmentation  and  develop  a  procedure  to  preserve  orthonormality  between  blocks.  This 
involves  adapting  an  developing  windowing  systems  for  time- frequency  analysis.  We  then  demonstrate 
the  connection  of  these  techniques  with  Gabor  and  wavelet  analysis. 

We  also  gave  a  technique  for  multi-rate  sub-Nyquist  sampling.  This  then  allows  for  processing  of 
wider  bandwidth  signals  at  smaller  bandwidth  rates.  These  techniques  are  based  on  our  work  on  multi¬ 
channel  deconvolution.  The  “tricks”  for  the  multi-rate  theory  come  from  clever  uses  of  number  theory 
and  complex  analysis. 

We  also  developed  a  program  of  study  connecting  sampling  theory  with  the  geometry  of  the  signal 
and  its  domain.  It  is  relatively  easy  to  demonstrate  this  connection  in  Euclidean  spaces,  but  one  quickly 
gets  into  open  problems  when  the  underlying  space  is  not  Euclidean.  In  Euclidean  space,  the  minimal 
sampling  rate  for  Paley-Wiener  functions  on  Mrf,  the  Nyquist  rate,  is  a  function  of  the  band-width.  No 
such  rate  has  yet  been  determined  for  hyperbolic  or  spherical  spaces.  We  look  to  develop  a  structure  for 
the  tiling  of  frequency  spaces  in  both  Euclidean  and  non-Euclidean  domains.  In  particular,  we  establish 
Nyquist  tiles  and  sampling  groups  in  Euclidean  geometry,  and  discuss  the  extension  of  these  concepts  to 
hyperbolic  and  spherical  geometry  and  general  orientable  surfaces. 

We  close  by  discussing  our  work  in  the  analysis  of  point  processes.  We  have  developed  a  collection  of 
algorithms  that  analyze  periodic  phenomena  generated  by  either  a  single  or  multiple  generators.  These 
work  even  when  the  data  is  extremely  sparse  and  noisy.  The  algorithms  use  number  theory  in  novel  ways 
to  extract  the  underlying  period(s)  by  modifying  the  Euclidean  algorithm  to  determine  the  period  from 
a  sparse  set  of  noisy  measurements.  We  divide  our  analysis  into  two  cases  periodic  processes  created 
by  a  single  source,  and  those  processes  created  by  several  sources.  We  wish  to  extract  the  fundamental 
period  of  the  generators,  and,  in  the  second  case,  to  deinterleave  the  processes. 
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1  Original  Statement  of  Objectives 

We  have  investigated  and  propose  to  continue  investigating  new  techniques  in  the  time- frequency  analysis 
of  signals.  The  new  techniques  are  developed  to  deal  with  classes  of  signals  for  which  the  known  techniques 
have  limitations,  e.g.,  adaptive  frequency  band  (AFB)  and  ultra- wide  band  (UWB)  signals.  These  signals 
are  of  interest  to  Air  Force  communications  systems,  but  the  theory  and  techniques  we  have  developed 
and  propose  to  continue  developing  fit  in  the  context  of  expanding  the  general  theory.  Our  techniques 
involve  signal  segmentation,  projection,  expansion  in  bases,  transform  analysis  and  multi-rate  methods. 

Adaptive  frequency  band  (AFB)  and  ultra-wide  band  (UWB)  systems  require  either  rapidly  chang¬ 
ing  or  very  high  sampling  rates.  Conventional  analog-to-digital  devices  have  non- adaptive  and  limited 
dynamic  range.  We  investigate  AFB  and  UWB  sampling  via  a  basis  projection  method.  The  method 
decomposes  the  signal  into  a  basis  over  time  segments  via  a  continuous-time  inner  product  operation 
and  then  samples  the  basis  coefficients  in  parallel.  The  signal  may  then  be  reconstructed  from  the  basis 
coefficients  to  recover  the  signal  in  the  time  domain.  We  develop  the  procedure  of  this  method,  analyze 
various  methods  for  signal  segmentation  and  develop  a  procedure  to  preserve  orthonormality  between 
blocks.  We  then  demonstrate  the  connection  of  these  techniques  with  Gabor  and  wavelet  analysis.  Items 
for  the  basis  projection  method  include: 

•  Develop  the  projection  method  for  a  variety  of  bases,  tailored  to  different  classes  of  input  signals. 
Target  signals  of  interest  to  Air  Force  communications  systems.  Implement  numerical  models  in 
MATLAB.  Propose  projection  systems  for  chip  design. 

•  Implement  Almost  Orthogonality  Theory  (Cotlar,  Knapp,  Stein)  to  develop  B-spline  windowing 
systems  that  are  relatively  easy  to  compute  and  almost  orthogonal  between  blocks. 

We  also  give  a  technique  for  multi-rate  sub-Nyquist  sampling.  This  then  allows  for  processing  of  wider 
bandwidth  signals  at  smaller  bandwidth  rates.  These  techniques  are  based  on  our  work  on  multi-channel 
deconvolution.  We  give  an  application  of  these  techniques  to  multi-rate  signal  streaming.  In  particular, 
we  will  investigate: 

•  Develop  a  systematic  approach  to  multi-rate  theory  in  order  to  develop  computable  systems. 

•  Construct  multi-rate  systems  for  signals  of  interest  to  Air  Force  communications.  Implement  nu¬ 
merical  models  in  MATLAB. 

We  close  by  placing  the  theories  in  the  context  of  standard  analysis  techniques.  We  also  look  into 
underlying  principles,  in  particular  for  sampling.  We  propose  developing  a  unified  approach  to  sampling 
via  “tiling  groups.”  This  allows  us  to  extend  sampling  into  non-Euclidean  geometries.  In  particular,  we 
outline  how  to  develop  sampling  for  hyperbolic  and  spherical  geometries.  Items  include: 

•  Develop  theory  of  “tiling  groups.” 

•  Use  the  theory  of  “tiling  groups”  to  develop  sampling  formulae  for  hyperbolic  and  spherical  geome¬ 
tries. 


2  Projects  Developed  During  the  Grant 

We  developed  new  techniques  and  theories  in  the  analysis  of  signals  over  the  duration  of  the  grant.  These 
techniques  and  theories  were  extensions  of  known  techniques  -  sampling,  Fourier,  Gabor  and  wavelet 
analysis,  and  new  approaches  to  analysis  -  using  combinations  of  analysis,  geometry,  group  theory,  and 
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number  theory.  Every  target  item  in  the  original  statement  of  objectives  was  achieved,  and 
several  new  areas  of  research  were  open  up.  There  were  four  main  areas  of  study. 

1. )  Techniques  were  developed  to  deal  with  classes  of  signals  for  which  the  known  techniques  have  limita¬ 

tions,  e.g.,  adaptive  frequency  band  (AFB)  and  ultra- wide  band  (UWB)  signals.  These  signals  are  of 
interest  to  Air  Force  communications  systems,  but  the  theory  and  techniques  we  have  developed  and 
propose  to  continue  developing  fit  in  the  context  of  expanding  the  general  theory.  Our  techniques 
involve  signal  segmentation,  projection,  expansion  in  bases,  transform  analysis,  and  multi-rate  meth¬ 
ods.  Adaptive  frequency  band  (AFB)  and  ultra-wide  band  (UWB)  systems  require  either  rapidly 
changing  or  very  high  sampling  rates.  Conventional  analog-to-digital  devices  have  non-adaptive  and 
limited  dynamic  range.  We  investigate  AFB  and  UWB  sampling  via  a  basis  projection  method.  We 
refer  to  the  technique  as  The  Projection  Method.  The  method  decomposes  the  signal  into  a  basis  over 
time  segments  via  a  continuous-time  inner  product  operation  and  then  samples  the  basis  coefficients 
in  parallel.  The  signal  may  then  be  reconstructed  from  the  basis  coefficients  to  recover  the  signal 
in  the  time  domain.  We  develop  the  procedure  of  this  method,  analyze  various  methods  for  signal 
segmentation  and  develop  a  procedure  to  preserve  orthonormality  between  blocks.  This  involves 
adapting  an  developing  windowing  systems  for  time- frequency  analysis.  We  then  demonstrate  the 
connection  of  these  techniques  with  Gabor  and  wavelet  analysis.  This  work  is  described  in  Section  4. 

2. )  We  also  gave  a  technique  for  multi-rate  sub-Nyquist  sampling.  This  then  allows  for  processing  of 

wider  bandwidth  signals  at  smaller  bandwidth  rates.  These  techniques  are  based  on  our  work  on 
multi-channel  deconvolution.  The  “tricks”  for  the  multi-rate  theory  come  from  clever  uses  of  number 
theory  and  complex  analysis.  This  work  is  described  in  Section  5. 

3. )  We  also  developed  a  program  of  study  connecting  sampling  theory  with  the  geometry  of  the  signal 

and  its  domain.  It  is  relatively  easy  to  demonstrate  this  connection  in  Euclidean  spaces,  but  one 
quickly  gets  into  open  problems  when  the  underlying  space  is  not  Euclidean.  In  Euclidean  space, 
the  minimal  sampling  rate  for  Paley-Wiener  functions  on  Mrf,  the  Nyquist  rate,  is  a  function  of  the 
band-width.  No  such  rate  has  yet  been  determined  for  hyperbolic  or  spherical  spaces.  We  look  to 
develop  a  structure  for  the  tiling  of  frequency  spaces  in  both  Euclidean  and  non-Euclidean  domains. 
In  particular,  we  establish  Nyquist  tiles  and  sampling  groups  in  Euclidean  geometry,  and  discuss  the 
extension  of  these  concepts  to  hyperbolic  and  spherical  geometry  and  general  orientable  surfaces. 
This  work  is  described  in  Section  6. 

4. )  We  close  by  discussing  our  work  in  the  analysis  of  point  processes.  This  work,  although  not  in  the 

original  proposal,  is  very  relevant  to  AFOSR  interests.  It  was  developed  over  the  course  of  the  grant. 

We  have  developed  a  collection  of  algorithms  that  analyze  periodic  phenomena  generated  by  either 
a  single  or  multiple  generators.  These  work  even  when  the  data  is  extremely  sparse  and  noisy.  The 
algorithms  use  number  theory  in  novel  ways  to  extract  the  underlying  period(s)  by  modifying  the 
Euclidean  algorithm  to  determine  the  period  from  a  sparse  set  of  noisy  measurements.  We  divide 
our  analysis  into  two  cases  periodic  processes  created  by  a  single  source,  and  those  processes  created 
by  several  sources.  We  wish  to  extract  the  fundamental  period  of  the  generators,  and,  in  the  second 
case,  to  deinterleave  the  processes. 

We  first  present  very  efficient  algorithm  for  extracting  the  fundamental  period  from  a  set  of  sparse 
and  noisy  observations  of  a  single  source  periodic  process.  The  procedure  is  computationally  straight¬ 
forward,  stable  with  respect  to  noise  and  converges  quickly.  Its  use  is  justified  by  a  theorem  which 
shows  that  for  a  set  of  randomly  chosen  positive  integers,  the  probability  that  they  do  not  all  share  a 
common  prime  factor  approaches  one  quickly  as  the  cardinality  of  the  set  increases.  The  proof  of  this 
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theorem  rests  on  a  probabilistic  interpretation  of  the  Riemann  zeta  function.  We  then  build  upon 
this  procedure  to  deinterleave  and  then  analyze  data  from  multiple  periodic  processes.  This  relies 
both  on  the  probabilistic  interpretation  of  the  Riemann  zeta  function,  the  equidistribution  theorem 
of  Weyl,  and  Wieners  periodogram.  This  work  is  described  in  Section  7. 


3  Introduction:  Background  and  Notation 

All  functions  considered  in  this  report  are  absolutely  and  square  integrable  functions  (/  G  L1  n  L2), 
unless  noted  otherwise.  References  for  the  material  on  harmonic  analysis  include  Benedetto  [Ben97], 
Daubechies  [Dau92]  Dym  and  McKean  [DM72],  Grafakos  [Gral4],  Higgins  [Hig96],  Homander  [Hor90], 
Levin  [Lcv96],  Papoulis  [Pap62,  Pap77],  Young  [YouSO]  and  Zayed  [Zay93].  Let  /  be  aperiodic,  integrable 
function  on  M,  with  period  2<L,  i.e. ,  /  G  L1(T2$).  The  Fourier  coefficients  of  /,  f[n],  are  defined  by 
f[n\  =  9^  /_$  f{t)  exp(— wrnt/<E>)  dt .  If  {/[n]}  is  absolutely  summable  ({/[n]}  G  l1),  then  the  Fourier 
series  of  /  is  f(t)  =  f[n\exp(iTrnt/&) .  Let  /  be  a  function  in  Ll.  The  Fourier  transform  of  / 

is  defined  as  /( co)  =  fR  f(t)e~2mtuldt  for  t  G  M  (time),  u>  G  M  (frequency).  The  inversion  formula, 
for  /  G  L1(M),  is  f(t)  =  (/)  (t)  =  Jg  f(tv)e2mujtduj.  Formally,  we  can  think  of  the  transform  and  the 
coefficient  integral  as  analysis,  and  the  inverse  transform  and  series  as  synthesis.  The  choice  to  have  27t 
in  the  exponent  simplifies  certain  expressions,  e.g.,  for  f,g  G  L1  n  L2(M),  /, g  G  L1  n  L2(M),  we  have 
the  Parseval  -  Plancherel  equations  -  ||/||l2(k)  =  II/IIl2(m)  and  (f,d)  =  (fid)  ■  We  extend  the  transform 
from  L 1  n  L2  to  L2  via  a  density  argument. 

Classical  sampling  theory  applies  to  square  integrable  band-limited  functions.  A  function  that  is  both 
H  bandlimited  and  L 2  has  several  smoothness  and  growth  properties  given  in  the  Paley- Weiner  Theorem 
(see,  e.g.,  [DM72]).  We  define  PWq  =  {/:/,/  G  L2,supp(/)  C  [— H,D]}.  The  Whittaker-Kotel’nikov- 
Shannon  (W-K-S)  Sampling  Theorem  applies  to  functions  in  PWq. 


Theorem  1  (W-K-S  Sampling  Theorem)  Let  f  G  PWq,  sine x{t)  =  sm^ ^  and  dnT(t)  =  5(t  —  nT). 
a.)  If  T  <  1/2 Q,  then  for  all  iGl, 


f(t)  =  Tj2f(nT) 

nEZ 


sin(^(t  —  nT)) 
7r  (t  —  nT) 


=  T 


LneZ 


dnT 


f 


*  sinc(t) 
T 


b.)  If  T  <  1/2 Q  and  f(nT)  =  0  for  all  n  G  Z,  then  f  =  0. 


4  The  Projection  Method 

Adaptive  frequency  band  (AFB)  and  ultra-wide  band  (UWB)  systems  are  of  interest  to  the  ASFOSR 
signal  processing  community.  They  require  either  rapidly  changing  or  very  high  sampling  rates,  and  this 
severely  stresses  classical  sampling  approaches.  At  UWB  rates,  conventional  analog-to-digital  devices 
have  limited  dynamic  range  and  exhibit  undesired  nonlinear  effects  such  as  timing  jitter.  Increased 
sampling  speed  leads  to  less  accurate  devices  that  have  lower  precision  in  numerical  representation.  This 
motivates  alternative  sampling  schemes  that  use  mixed-signal  approaches,  coupling  analog  processing 
with  parallel  sampling,  to  provide  improved  sampling  accuracy  and  parallel  data  streams  amenable  to 
lower  speed  (parallel)  digital  computation.  Developing  a  sampling  theory  for  adaptive  frequency  band 
and  ultra-wide  band  systems  is  the  motivation  for  the  methods  in  this  section  of  the  report.  Two  of  the 
key  items  needed  for  this  approach  are  quick  and  accurate  computations  of  Fourier  coefficients,  which 
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are  computed  in  hardware  and  for  which  we  have  developed  a  fast  scheme  ([Casl2,  CC15,  CaslS]),  and 
effective  adaptive  windowing  systems  ([Cast  ,  ,  ]).  Sampling  schemes  that  work  with  inverse 

operators  and  multi-rate  schemes  were  developed  in  [CW94,  Cas98,  CS99a,  CW01,  Cas07]. 

We  investigate  AFB  and  UWB  sampling  via  a  basis  projection  method.  The  method  was  introduced 
as  a  means  of  UWB  parallel  sampling  and  applied  to  UWB  communications  systems.  The  method  first 
decomposes  the  signal  into  a  basis  over  time  segments  via  a  continuous-time  inner  product  operation 
and  then  samples  the  basis  coefficients  in  parallel.  The  signal  may  then  be  reconstructed  from  the  basis 
coefficients  to  recover  time  domain  samples,  or  further  processing  may  be  carried  out  in  the  new  domain. 
We  address  several  issues  associated  with  the  basis-expansion  and  sampling  procedure,  including  choice 
of  basis  and  segmentation  of  the  signal.  We  develop  a  mathematical  model  of  the  procedure,  using  both 
standard  (sine,  cosine)  basis  elements  and  general  basis  elements,  and  give  this  representation  in  both 
the  time  and  frequency  domains.  First,  we  consider  ideal  rectangular  windowing  of  each  segment,  and 
then  develop  a  theory  for  overlapping  windowed  segments.  Splines  are  used  to  analyze  the  degree  of 
smoothness  in  time  and  corresponding  decay  in  frequency  that  occurs  with  overlapping  windows.  We 
employ  the  theory  of  lapped  orthogonal  transforms  to  preserve  the  orthogonality  of  basis  elements  in 
the  overlapping  regions.  We  compare  the  method  with  traditional  sampling  and  close  by  applying  the 
procedure  to  binary  signals.  Proofs  of  all  results  and  extensive  error  analysis  of  the  projection  method 
can  be  found  in  Casey  and  Sadler  [Casl2]. 

The  theory  of  splines  and  some  techniques  from  ordinary  differential  equations  give  us  the  tools  to 
cut  up  the  time  domain  into  perfectly  aligned  segments  so  that  there  is  no  loss  of  information.  We  want 
the  systems  to  be  smooth,  so  as  to  provide  control  over  decay  in  frequency,  have  variable  cut-off  functions 
for  ffexibility  in  design,  and  adaptive,  so  as  to  adjust  accordingly  to  changes  in  frequency.  We  also  want 
to  develop  our  systems  so  that  the  orthogonality  of  bases  in  adjacent  and  possible  overlapping  blocks  is 
preserved. 


Definition  1  (ON  Window  System)  Let  0  <  r  <C  T.  An  ON  Window  System  is  a  set  of  functions 
(Wfc(f)}  such  that  for  all 

(i.)  supp(Wfc(t))  C  [kT  —  r,(k  +  1  )T  +  r]  , 

(ii.)  W k{t)  =  1  for  t  6  [kT  +  r,  (, k  +  1  )T  —  r\  , 

(in.)  W k((kT  +  T/2)  -  t)  =  W k(t  -  (kT  +  T/2))  for  t  G  [0,  T/2  +  r] , 

(iv.)  ^[Wfc(f)]2  =  l, 

(v.)  {[Wfc]°[n]}  is  absolutely  summable,  i.e.  {[Wfc]°[n]}  G  l1  ■ 


Conditions  (i.)  and  (ii.)  are  partition  properties,  in  that  they  give  an  exact  snapshot  of  the  input 
function  /  on  [kT  +  r,  (k  +  1  )T  —  r],  with  smooth  roll-off  at  the  edges.  Condition  (in.)  is  needed  to 
preserve  orthogonality  between  adjacent  blocks.  Condition  (iv.)  simplifies  computations,  and  condition 
(v.)  is  needed  for  the  computation  of  Fourier  coefficients. 

Let  {<Pj(t)}  be  an  orthonormal  basis  for  L2[— T/2, T/2],  Define 


ifj(t) 


0 

<pj(t) 

-<Pj(-T-t) 

<Pj(T-t) 


|f|  >  T/2  +  r 
\t\  <  T/2 

-T/2  —  r  <  t  <  -T/2 
T/2  <  t  <  T/2  +  r  . 


We  refer  to  {<fij(t)}  as  the  folded  basis  elements.  Let  Ta  be  the  translation  operator,  i.e.,  Ta[f](t)  = 
f(t  —  a)  .  All  of  the  basis  elements  are  presented  in  the  interval  [T /2  —  r,  T/2  +  r]  centered  at  the  origin. 
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Therefore,  the  operator  T\(k)T+T /2\  will  place  the  basis  in  the  interval  [( k)T  —  r,(k  +  1  )T  +  r].  In  the 
following,  W /(f)  is  the  window  centered  at  the  origin,  and  p>j  is  a  basis  element  in  that  window.  Let 

=  {^[(fc)T+T/2]  [W lPj\(t)}  . 

We  denote  this  as  {’Lfcj}  =  {Wfc</?)(t)} 

Theorem  2  (The  Orthogonality  of  Overlapping  Blocks)  The  collection 
{'I'fcj}  =  {Wfc^j(t)}  forms  an  orthonormal  basis  for  L2(M). 

Given  characteristics  of  the  class  of  input  signals,  the  choice  of  basis  functions  used  in  the  projection 
method  can  be  tailored  to  optimal  representation  of  the  signal  or  a  desired  characteristic  in  the  signal. 


Theorem  3  (Projection  Formula  for  ON  Windowing)  Let  {Wfc(f)}  be  an  ON  Window  System, 
and  let  {’Lfcj}  be  an  orthonormal  basis  that  preserves  orthogonality  between  adjacent  windows.  Let 
f  G  PWq  and  N  =  N(T,Ll)  be  such  that  {f,^k,n}  =  0  for  all  n  >  N  and  all  k.  Then,  f(t)  f-p(t), 
where 


bit)  =  ^2 


r  N 

E 


(/,  *k,n)*k,n(t) 


k£Z  Ln=-N 


Given  the  flexibility  of  our  windowing  system,  we  can  also  formulate  an  adaptive  projection  system 
for  the  ON  Window  Systems. 

Theorem  4  (Adaptive  Projection  Formula  for  ON  Windowing)  Let  f,  f  G  L2(M)  and  f  have 
a  variable  but  bounded  band-limit.  Ll(t).  Let  r(t)  be  an  adaptive  block  of  time.  Let  {Wfc(f)}  be  an  ON 
Window  System  with  window  size  r(t)  +  2 r  on  the  kth  block,  and  let  {'Lfcj}  be  an  orthonormal  basis 
that  preserves  orthogonality  between  adjacent  windows.  Given  r(t),  let  Ll(t)  =  max{fl(f)  :  t  G  r(f)}.  Let 
N(t)  =  N(r(t),  Ll(t))  be  such  that  {f,^k,n)  =  0.  Then,  f(t)  ss  fp{t),  where 

r  N(t)  n 

M*)=E  E  (f^k,n)^k,n(t)  ■ 

&EZ  n=—N(t) 


4.1  Remarks  on  Applications 

Despite  extensive  advances  in  integrated  circuit  design  and  fabrication  processes,  wideband  problems 
continue  to  hit  barriers  in  sample  and  hold  architectures  and  analog-to-digital  conversion  (ADC).  ADC 
signal-to-noise  and  distortion  ratio,  the  effective  number  of  resolution  bits,  declines  with  sampling  rate 
due  to  timing  jitter,  circuit  imperfections,  and  electronic  noise.  ADC  performance  (speed  and  total 
integrated  noise)  can  be  improved  to  some  extent,  e.g.,  by  cooling  and  therefore  lowering  the  system 
temperature.  However,  the  energy  cost  may  be  significant,  and  this  presents  a  major  hurdle  for  imple¬ 
mentation  in  miniaturized  devices.  Digital  circuitry  has  provided  dramatically  enhanced  digital  signal 
processing  operation  speeds,  but  there  has  not  been  a  corresponding  dramatic  energy  capacity  increase 
in  batteries  to  operate  these  circuits;  there  is  no  Moore’s  Law  for  batteries  or  ADCs. 

A  growing  number  of  applications  face  this  challenge,  such  as  miniature  and  hand-held  devices  for 
communications,  robotics,  and  micro  aerial  vehicles.  Very  wideband  sensor  bandwidths  are  desired  for 
dynamic  spectrum  access  and  cognitive  radio,  radar,  and  ultra-wideband  systems.  Multi-channel  and 
multi-sensor  systems,  array  processing  and  beamforming,  multi-spectral  imaging,  and  vision  systems 
compound  the  issue.  All  of  these  rely  on  analog  sensing  and  a  digital  interface,  perhaps  with  feedback. 
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This  motivates  mixed-signal  circuit  designs  that  tightly  couple  the  analog  and  digital  portions,  and  operate 
with  parallel  reduced  bandwidth  paths  to  relax  ADC  requirements.  The  goal  of  such  wideband  integrated 
circuit  designs  is  to  achieve  good  tradeoffs  in  dynamic  range,  bandwidth,  and  parallelization,  while 
maintaining  low  energy  consumption.  This  requires  a  careful  balance  of  analog  and  digital  functionality. 

From  a  signal  processing  perspective,  we  have  approached  this  problem  by  implementing  an  appro¬ 
priate  signal  decomposition  in  the  analog  portion  that  provides  parallel  outputs  for  integrated  digital 
conversion  and  processing.  This  naturally  leads  to  an  architecture  with  windowed  time  segmentation 
and  parallel  analog  basis  expansion.  In  this  part  of  the  report,  we  viewed  this  from  the  sampling  theory 
perspective,  including  segmentation  and  window  design,  achieving  orthogonality  between  segments,  basis 
expansion  and  choice  of  basis,  signal  filtering,  and  reconstruction.  The  approach  we  have  developed  is  tai¬ 
lored  toward  strong  connections  to  circuit  design  considerations  and  applications.  We  look  to  generalize 
the  windowing  systems  to  obtain  computationally  efficient  approaches  to  time-frequency  analysis. 


5  Multi-rate  Sampling 


A  key  step  in  solving  the  multichannel  deconvolution  problem  involves  solving  an  interpolation  problem, 
reconstructing  functions  (the  deconvolvers)  in  a  space  of  restricted  growth  (S')  from  discrete  data  (their 
values  on  the  zero  sets  of  the  convolvers).  This  gives  solutions  to  the  Bezout  equation.  This  develop¬ 
ment  utilizes  the  zero  sets  of  the  /!,;  as  different  sampling  rates.  This  connection  between  sampling  and 
multichannel  deconvolution  naturally  leads  to  sampling  schemes  on  properly  chosen  non-commensurate 
lattices.  Let  {ri}rfL1  be  a  set  such  that  ri/rj  is  poorly  approximated,  by  rationals  for  i  /  j.  and  let 

R  =  YT=  iri-  Let 

A=u{t~}  urn 

i=1  l  )  fceN 


be  a  sampling  grid,  made  up  of  a  union  of  sampling  grids  with  non-commensurate  generators  {r*}.  Given 
a  R  bandlimited  function  /,  we  can  reconstruct  /  from  samples  of  /  taken  on  A.  We  refer  to  this  as  the 
Multi-Rate  Sampling  Problem  or  MRSP.  Our  work  on  this  can  be  found  in  [Casl6,  CSOO,  CS99a,  CS99b, 
CW01], 

Let  t  E  M,  and  let  a  be  an  irrational  that  is  poorly  approximated  by  rationals,  e.g.,  \/2,  or  (l  +  \/5)/2. 
Let  Hi  (t)  =  X[_1;1]  (t)  ,  p2  (t)  =  A[_Q  Q,](t)  model  the  impulse  response  of  the  channels  of  a  two-channel 
system.  Then 


Mi(C) 


sin(27r£) 

vrC 


/L>(C) 


sin(27ra£) 

vrC 


Since  a  is  poorly  approximated  by  rationals,  {pi}  is  strongly  coprime.  Now  let 


A,  = 


,  A-2  = 


for  n,m  E  N,  and  let  A  =  Ai  U  A2  .  Note  that  the  information  contained  in  the  original  signal  is 
reconstructed  by  creating  deconvolvers  defined  initially  on  A  U  {0}.  Order  the  elements  of  A,  denoting 
this  as  A  =  {A^}.  Walnut  (see  [Cash  ])  has  shown  the  following. 

Theorem  5  Let  a  be  an  irrational  that  is  poorly  approximated  by  rationals,  and  let  f  be  a  (1  +  ^-band- 
limited  function.  Then  f  is  uniquely  determined  by  {/(A*,)}  U  {/( 0),  /'( 0)}  . 


Note,  for  a  (1  +  cc)-band-limited  function,  the  Nyquist  rate  is  1  / (2(1  +  a)).  However,  our  individual 
sampling  rates  are  1/2  and  l/(2ce).  Both  these  rates  are  below  Nyquist.  The  reconstruction  of  /  from  this 
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lattice  is  achieved  by  using  complex  interpolation  theory.  These  techniques  go  back  to  basic  Lagrange 
interpolation,  and  were  developed  for  entire  functions  by  various  mathematicians,  most  notably  B.  Ya. 
Levin  [  ].  Let  G(z)  =  sin(27r,z)  sin(27raz)  .  Then,  G(z)  is  an  entire  function,  which  is  almost  periodic 

on  M,  and  has  simple  zeros  on  A,  and  a  double  zero  at  {0}. 

Proposition  1  Let  Aj,  G  A  and  let 

H(z)= _ m _ 

,U  G'(\j){z-\j)  ' 

Then  Hj( Xk)  =  6jk  . 

At  2  =  0,  we  have  to  construct  interpolating  functions  K\ ,  K2  so  that 

Ah(0)  =  1 ,  K[( 0)  =  0  ,  K2( 0)  =  0  ,  K'2( 0)  =  1 . 

Using  the  Taylor  expansion  of  G,  we  derive  the  following. 


Proposition  2 

Ki(z)  = 


G(z) 


G(z) 


,  K2(z)  = 


G{z) 


G{z ) 


(G"  (fd) /2\)(z2)  (47r2a)(z2)  ’  (G"(X)/2\)(z)  (47r2a)(^) 

Combining  these  two  propositions  gives  us  the  reconstruction  formula. 

Theorem  6  Let  a  be  an  irrational  that  is  poorly  approximated  by  rationals,  and  let  f  be  a  (1  +  ^-band- 
limited  function.  Let 


Ai  = 


(±k 

l“2“ 


,  A2  = 


±k\ 

2 a  J  ’ 


for  k  £  N,  and  let  A  =  {\k}  =  Ai  U  A2  .  Then  f  is  uniquely  determined  by  {f(Xk)}  U  {/(0),  /7(0)}  . 
Moreover,  f  can  be  approximated  from  its  values  on  A  U  {0}  by  the  formula 

+;(0)  cw +m 


where 


(47r2a)(t2)  '  J  y  y  (47 r2a)(t) 

G(t)  =  sin(27rt)  •  sin(2a7rt) . 


Convergence  is  uniform  convergence  on  compact  subsets,  and  is  shown  in  a  manner  similar  to  the 
results  for  deconvolution.  The  proof  is  lengthy,  and  appears  in  [  asl6].  We  need  to  point  out  two 
items.  First,  the  sampling  grid  is  rigid.  Perturbation  of  the  grid  results  in  a  loss  of  information.  Sec¬ 
ond,  that  because  sampling  points  in  A  can  get  arbitrarily  close  together  (inf{|Am  —  An|}  =  0,  and  so 
inf  irifcYj  I  sin(27rrfcA)|  :  A  E  A^}  =  0),  the  set  of  interpolating  functions  can  not  form  a  Riesz  basis  and 
the  interpolating  formula  can  not  converge  in  norm  (see  [  ]).  In  fact,  the  interpolating  functions  do 

form  a  Bessel  sequence,  but  do  not  form  a  frame  and  therefore  do  not  form  a  Riesz  basis.  The  problems 
occur  at  points  where  the  sampling  points  get  close  together.  The  interpolating  function  follows  the  orig¬ 
inal  function  along  exactly,  except  for  a  very  subtle  “ripple”  at  those  points  where  the  sampling  points 
get  close  together.  The  following  stabilizes  the  construction.  Let  d  be  given,  0  <  <5  <  Let 

A<5  =  {A  E  A  :  dist(A,  A  \  {A})  <  5}  . 
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Elements  in  As  occur  in  pairs,  with  each  pair  containing  one  element  from  Ai  =  {^r}  and  the  other  from 
A2  =  {  }  for  k  G  N.  Let  {Afc}  =  A  =  Ai  U  A2  =  As  U  A  /  ■  For 


G(t)  =  sin(27rf)  •  sin(27rat),  let  H\k  = 


G'{Xk)(t-Xk)  ' 


We  can  show  that 


G{t) 

G'{\k){t  -  Afc) 


is  dense  in  Li,  .  ,  and  that 
(!+«) 


f  G(t )  \ 

{G>(\k)(t-\k)jXk£As, 

is  a  Bessel  sequence.  Also,  we  can  show  that 


G(t)  G(t) 

47T at  ’  47rat2 


6  Closure  (  Span 


G(t) 

G’(Xk)(t-Xk) 


The  points  in  As  can  be  treated  as  double  points.  For  { Ai ,  A2}  C  As, 

[f(Xi)HXl(t)  +  f(X2)HX2(t)] 


=  /(A1)ifXl(t)  +  /'«Al)[7— — — ^4 - 

1  (  1 -  Amt  ~  a2>  J 

+  7£(Ai  —  A2)2  , 

where  £\i  — >  Ai  and  1Z(X\  —  A2)2  — >  0  as  <5  — >  0.  Finally,  we  can  show  that 

f  G(t)  \ 

l  {G”(Xk)/2\)(t  -  Afc)2  J  \keAs 

is  a  Bessel  sequence. 

The  result  generalizes.  We  can  create  sampling  sets  on  £  lattices  using  a  set  of  numbers  {rt}j=1  such 
that  (r'i/rj)  is  poorly  approximated  by  rationals  for  i  ^  j.  Again,  if  is  a  set  of  primes, 

{1,  \/pT,  y/PlP2,  ■  ■  ■,  y/PlP2  ■■■pe-l} 

is  a  set  of  numbers  whose  ratios  are  poorly  approximated  by  rationals.  Let  Afc  =  for  n£N,  and 

let 


We  reconstruct  on  A  U  {0},  letting 


IJ  Afc  =  A  =  {A/J  . 


c(u = n  sin(27r?’fcz) 


and  letting 


H  (z)  =  _ M _ 

mU  G'(Xm)(z~Xm) 
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Then  Hm(Xn )  =  5mn  .  The  interpolating  functions  at  the  origin  are  a  linear  combination  of 

G(z)  , 


9 j  = 


(*»') 


chosen  so  that 


Klj  1)(0)  =  6kj,k,j  =  l,...,e. 

As  before,  using  Taylor  series,  we  get 


2n+27rn+2  („_fc  +  2)!  _ 

3!  (n-  A;-j  +  2)!A/Pi:P2‘“:Pri 

71—1 

(1  +  J>7)in-fc-j  +  2 

+0(tn-k-j+4)  _ 


We  want 

k 

Kk[t)  =  ^2  CiHi(t )  so  that  ^_1)(0)  =  j  =  1,  •  •  ■ ,  i  ■ 

1=1 

Solving  these  relationships  gives 

k 

Kk(t)  = 

i= 1 

E,rtot|/21(Sr(i  + 

■  (fc  -  1)!  • 


6  Sampling  and  Non-Euclidean  Geometry 

The  study  of  non-Euclidean  geometries  is  becoming  increasingly  relevant  in  signal  processing.  There 
are  numerous  motivations  for  extending  signal  processing,  and  in  particular,  sampling  theory,  to  non- 
Euclidean  spaces,  and  in  particular,  hyperbolic  and  spherical  geometries.  Hyperbolic  space  and  its  im¬ 
portance  in  Electrical  Impedance  Tomography  (EIT)  [  CT96,  BerOl]  and  Network  Tomography  [:  1GB06] 
has  been  mentioned  in  several  papers  of  Berenstein  et.  al.  and  some  methods  developed  in  papers  of 
Kuchment,  e.g.,  [KucOG].  Irregular  sampling  of  band-limited  functions  by  iteration  in  hyperbolic  space 
is  possible,  as  shown  by  Feichtinger  and  Pesenson  [  ,  ’P05]  and  Christensen  and  Olafsson  [  013]. 

Applications  where  data  are  defined  inherently  on  the  sphere  are  found  in  computer  graphics,  planetary 
science,  geophysics,  quantum  chemistry,  and  astrophysics  [DH94,  MM  ].  In  many  of  these  applications, 
a  harmonic  analysis  of  the  data  is  insightful.  For  example,  spherical  harmonic  analysis  has  been  remark¬ 
ably  successful  in  cosmology,  leading  to  the  emergence  of  a  standard  cosmological  model  [CW98,  MW11]. 

Our  approach  is  to  connect  sampling  theory  with  the  geometry  of  the  signal  and  its  domain  [  ’  C 1 6] . 
It  is  relatively  easy  to  demonstrate  this  connection  in  Euclidean  spaces,  but  one  quickly  gets  into  open 
problems  when  the  underlying  space  is  not  Euclidean.  For  Paley-Wiener  functions  on  M  the  minimal 
sampling  rate,  the  Nyquist  Rate ,  is  a  function  of  the  band-width.  The  establishment  of  the  exact  Nyquist 
rate  in  non-Euclidean  spaces  is  an  open  problem.  This  rate  is  needed  to  develop  regular  sampling.  We 
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use  two  tools  to  work  on  the  problem  -  the  Beurling-Landau  density  and  Voronoi  cells.  Using  these  tools, 
we  establish  a  relation  in  Euclidean  domains,  connecting  Beurling-Landau  density  to  sampling  lattices 
and  hence  dual  lattice  groups,  and  then  use  these  dual  lattices  to  define  Voronoi  cells,  which  become  our 
tiles  in  frequency.  We  then  discuss  how  to  extend  this  to  non-Euclidean  geometries. 

Classical  sampling  theory  applies  to  square  integrable  band-limited  functions.  A  function  that  is  both 
bandlimited  and  L 2  has  several  smoothness  and  growth  properties  given  in  the  Paley- Weiner  Theorem 
(see,  e.g.,  [DM72]).  We  define  PWq  =  {/:/,/  G  L2,supp(/)  C  [— U,U]}.  The  Whittaker-Kotel’nikov- 
Shannon  (W-K-S)  Sampling  Theorem  applies  to  functions  in  PWq. 


Theorem  7  (W-K-S  Sampling  Theorem)  Let  f  G  PWfj,  sincT(i)  =  ^  and  dnT(t )  =  5(t  —  nT ). 

a.)  If  T  <  1/2U,  then  for  all  iGl, 


f(t)=Tj2f(nT) 

n£L 


sin  (^(t  —  nT )) 
7r  (t  —  nT) 


=  T 


1 ZSnT 


SiEZ 


/  )  *  sinc(f) 


b.)  If  T  <  1/2 Q  and  f(nT)  =  0  for  all  rGZ,  then  f  =  0. 

The  sphere  §2  is  compact,  and  its  study  requires  different  tools.  Fourier  analysis  on  S2  amounts 
to  the  decomposition  of  L2(§2)  into  minimal  subspaces  invariant  under  all  rotations  in  50(3) .  Band- 
limited  functions  on  the  sphere  are  spherical  polynomials.  Sampling  on  the  sphere  is  how  to  sample  a 
band-limited  function,  an  IVth  degree  spherical  polynomial,  at  a  finite  number  of  locations,  such  that 
all  of  the  information  content  of  the  continuous  function  is  captured.  Since  the  frequency  domain  of  a 
function  on  the  sphere  is  discrete,  the  spherical  harmonic  coefficients  describe  the  continuous  function 
exactly.  A  sampling  theorem  thus  describes  how  to  exactly  recover  the  spherical  harmonic  coefficients 
of  the  continuous  function  from  its  samples.  Developing  sampling  lattices  leads  to  questions  on  how  to 
efficiently  tile  the  sphere,  a  subject  in  its  own  right,  e.g.,  Driscoll  and  Healy  [DH94],  Keiner,  Kunis,  and 
Potts  [KKP07],  and  McEwen  and  Wiaux  [V 1 '  ]. 

The  Nyquist  rate  allows  us  to  develop  an  efficient  tiling  of  frequency  space.  A  tiling  or  a  tessellation 
of  a  flat  surface  is  the  covering  of  the  plane  or  region  in  the  plane  using  one  or  more  geometric  shapes, 
called  tiles,  with  no  overlaps  and  no  gaps.  This  generalizes  to  higher  dimensions.  We  look  to  develop 
Nyquist  tiles  and  sampling  groups  for  Euclidean,  hyperbolic,  and  spherical  spaces  and  general  analytic 
surfaces.  We  first  assume  that  signals  are  single  band  and  symmetric  in  frequency,  i.e.,  that  the  transform 
of  the  signal  can  be  contained  in  a  simply  connected  region  centered  at  the  origin.  Symmetry  can  be 
achieved  by  shifting,  and  multi-band  signals  can  be  addressed  by  the  techniques  in  this  report. 

Consider  sampling  in  M.  Given  /  G  PWq  choose  T  <  1/(20).  For  /  G  L1([0,T))  the  T -periodization  of 
/  is  (/t)°(£)  =  Yfn&z  f(t~nT) .  We  can  then  expand  (/r)°(t)  in  a  Fourier  series.  The  sequence  of  Fourier 

coefficients  of  this  T-periodic  function  are  given  by  {fT)°[n\  =  p  f  (— •  Then  if  the  value  of  the  Fourier 
series  of  fr°  agrees  with  the  value  of  /t°  at  the  origin,  we  get  T  Yhn&i  f(n^)  =  Snez  /(n/T)PSFl . 
Thus,  the  Poisson  Summation  Formula  allows  us  to  compute  the  Fourier  series  of  (/t)°  in  terms  of  the 
Fourier  transform  of  /  at  equally  spaced  points.  This  extends  to  the  Schwartz  class  of  distributions  as 
$ nT  =  y  zFnez  ^n/TPSF2  .  Then,  for  /  G  PWq,  if  we  sample  at  exactly  Nyquist 


m  = 


2U 


2>(^) 

LnEZ 


*  sinc(t) 

U) 
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if  and  only  if 


/H 


E 

AnGZ  L 


^2  nf2 


/  -x. 


Q,Q)(W)  • 


The  interval  [— f2, 11)  is  simply  connected  and  symmetric  to  the  origin.  It  is  spread  by  the  group  of 
translations  to  form  a  tiling  of  frequency  space  -  {[(A;  —  1)11,  (k  +  1)H)}.  We  refer  to  [— 11,11)  as  a 
sampling  interval.  Note,  sampling  intervals  are  “half  open,  half  closed,”  with  length  determined  by  the 
Nyquist  rate.  The  inverse  transform  of  the  characteristic  functions  of  the  tiles  are  sine  functions,  which 
form  an  orthonormal  (o.n.)  basis  for  PWq.  Sampling  is  expressed  in  terms  of  this  basis.  We  can  now 
define  the  following. 


Definition  2  (Nyquist  Tiles  for  /  £  PWq)  Let  f  be  a  non-trivial  function  in  PWq.  The  Nyquist  Tile 
NT(/)  for  f  is  the  sampling  interval  of  minimal  length  in  M  such  that  supp(/)  C  NT(/) .  A  Nyquist  Tiling 
for  f  is  the  set  of  translates  {NT(/)fc}fcg^  of  Nyquist  tiles  which  tile  M. 

The  Nyquist  tile  is  transported  by  a  group  of  motions  to  cover  the  transform  domain. 


Definition  3  (Sampling  Group  for  /  £  PWq)  Let  f  £  PWq  with  Nyquist  Tile  NT(/).  The  Sampling 
Group  G (/)  is  a  group  of  translations  such  that  NT (/)  tiles  M. 


We  now  extend  these  ideas  to  A-dimensional  Euclidean  space.  Let  T  >  0  and  let  fit)  be  a  function 
such  that  supp /  C  [0,T]fc.  The  T- periodization  of  /  is  [/]°(t)  =  Ylnezd  fit  —  nT)  .  We  can  expand 
a  T-periodic  function  [/]°(f)  in  a  Fourier  series.  Denote  the  lattice  A  =  TZd,  where  T  is  the  n  x  n 
matrix  with  T  on  the  main  diagonal  and  zeroes  elsewhere.  The  sequence  of  Fourier  coefficients  of  this 
periodic  function  on  the  lattice  A  =  TZrf  are  given  by  [/]°[n]  =  ^d/(-f)  •  We  have  +  nT)  = 

T?  Enezd  f{n/T)e2'Kin'tlT  .  Therefore,  Ene lA  f{nT)  =  ^  Enez^  f(n/T)FSFl .  We  can  write  the  Poisson 
summation  formula  for  an  arbitrary  lattice  by  a  change  of  coordinates.  Let  A  be  an  invertible  dxd  matrix, 
A  =  AZd,  and  A1-  =  (A7  )-1Zd  be  the  dual  lattice.  Then 

E  /(*  +  *)  =  E(/°AHA_lt  +  ">  =  Ew  E  /((A Tr1(n))e2"<ATU,<">*. 

AgA 


Note,  |detA|  =  vol(A).  This  last  expression  can  be  expressed  more  directly  as  EasA  /(^  +  A)  = 
vQj~  E/3eAx  /(/3)e2,r*^'t .  This  extends  again  to  the  Schwartz  class  of  distributions  as 

The  sampling  formula  again  follows  from  computations  and  an  application  of  (PSF2).  We  assume 
a  single  band  signal.  Let  A  be  a  regular  sampling  lattice  in  Md,and  let  A^  be  the  dual  lattice  in  Mrf. 
Then  A  has  generating  vectors  {ri,  T2, . . . ,  r^},  and  the  sampling  lattice  can  be  written  as  A  =  {A  : 
A  =  z\T\  +  Z2T 2  +  . . .  +  ZdTd}  for  (zi,  Z2,  ■  ■  ■ ,  Zd)  £  Zd.  Let  {fii,  Tto,  ■  ■  ■ ,  be  the  generating  vectors 
for  the  dual  lattice  A-1.  The  dual  sampling  lattice  can  be  written  as  A1-  =  {A^  :  A-1  =  zifl\  +  22^2  + 
. . .  +  ZdTtd}  for  (zi,  Z2,  ■  ■  ■  ,Zd)  £  Zd.  The  vectors  {fXi,  D2,  •  •  • ,  generate  a  parallelepiped.  We  want 
to  use  this  parallelepiped  to  create  a  tiling,  and  therefore  we  make  the  parallelepiped  “half  open,  half 
closed”  as  follows.  If  we  shift  the  parallelepiped  so  that  one  vertex  is  at  the  origin,  we  include  all  of 
the  boundaries  that  contain  the  origin,  and  exclude  the  other  boundaries.  We  denote  this  region  as  a 
sampling  parallelepiped  fl-p.  If  the  region  Lip  is  a  hyper-rectangle,  we  get  the  familiar  sampling  formula 

1  w-  „,ni  n^sin(^(f-mcui))  sinj^jt  -  ndwd)) 
wi""’  uid  7r(t  -  n±LOi)  ir{t  —  ndud) 


m  = 


vol(A) 


E  /( 
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If,  however,  the  sampling  parallelepiped  D-p  a  general  parallelepiped,  we  first  have  to  compute  the  inverse 
Fourier  transform  of  Xnp.  Let  S  be  the  generalized  sine  function  S  =  vq^  (XnP)v  .  Then,  the  sampling 

formula  (see  [Hig96])  becomes  f(t)  =  /(A)«S(t  ~  • 

Definition  4  (Nyquist  Tiles  for  /  G  PWnP)  Let 

TWfjp  =  {/  continuous  :  /  €  L2(Rd),  f  G  L2(Rd),  supp(/)  C  Dp}  , 

where  (Di,  D2, . . . ,  Lid}  be  the  generating  vectors  for  the  dual  lattice  A2-.  Let  f  be  a  non-trivial  function 
in  PWjjp-  The  Nyquist  Tile  NT (/)  for  f  is  the  sampling  parallelepiped  of  minimal  area  in  Rd  centered  at 
the  origin  such  that  supp(/)  C  NT (/).  A  Nyquist  Tiling  is  the  set  of  translates  {NT (f)k}k£Zd  of  Nyquist 
tiles  which  tile  Rd. 


Definition  5  (Sampling  Group  for  /  G  PWnp)  Let  f  G  PWnP  with  Nyquist  Tile  NT(/).  The  Sam¬ 
pling  Group  G  is  a  symmetry  group  of  translations  such  that  NT(/)  tiles  Rd. 


A  sequence  A  is  separated  or  uniformly  discrete  if  q  =  inffc(A k+i  —  A k)  >  0.  The  value  q  is  referred  to 
as  the  separation  constant  of  A.  With  a  separated  sequence  A  we  associate  a  distribution  function  nA(t) 
defined  such  that  for  a  <  b,  nA(b)  —  nA(a)  =  card(An  (a,  b\) ,  and  normalized  such  that  nA(0)  =  0.  There 
is  clearly  a  one-to-one  correspondence  between  A  and  nA.  A  discrete  set  A  is  a  set  of  sampling  for  PWq 
if  there  exists  a  constant  C  such  that  ||/|| 2  <  AfegA  |/(Afc)|2  for  every  /  G  PWq.  The  set  A  is  called  a 
set  of  interpolation  for  PWq  if  for  every  square  summable  sequence  (aA}AeAj  there  is  a  solution  /  G  PWq 
to  /(A)  =  a\,  A  G  A.  Clearly,  all  complete  interpolating  sequences  are  separated.  Landau  showed  that 
if  A  is  a  sampling  sequence  for  PWq,  then  there  exists  constants  A  and  B,  independent  of  a,  b  such  that 
nA(b)  —  nA(a)  >  [b  —  a)  —  Alog+(b  —  a)  —  B  . 


Definition  6  (Beurling-Landau  Densities)  The  Beurling-Landau  lower  D  and  upper  D+  densities 
are  given  by 


D  (A)  =  liminfJwooinfteK 


(nA(t  +  r))  -nA(t) 
r 


D+{  A)  =  lim  supr_>00supte]K 


(nA(t  +  r))  -nA(t) 
r 


The  densities  are  defined  similarly  in  higher  dimensions.  Specifically,  for  the  exact  and  stable  re¬ 
construction  of  a  band-limited  function  /  from  its  samples  {/(A*,)  :  A&  G  A},  it  is  sufficient  that  the 
Beurling-Landau  lower  density  satisfies  D~{ A)  >  1.  A  set  fails  to  be  a  sampling  set  if  D~( A)  <  1.  Con¬ 
versely,  if  /  is  uniquely  and  stably  determined  by  its  samples  on  A,  then  D~{ A)  >  1.  Note,  a  sampling  set 
for  which  the  reconstruction  is  stable  in  this  sense  is  called  a  (stable)  set  of  sampling.  This  terminology  is 
used  to  contrast  a  set  of  sampling  with  the  weaker  notion  of  a  set  of  uniqueness.  A  is  a  set  of  uniqueness 
for  PWq  if  / |A  =  0  implies  that  /  =  0.  Whereas  a  set  of  sampling  for  PWq  has  a  density  D~{ A)  >  1, 
there  are  sets  of  uniqueness  with  arbitrarily  small  density.  We  also  have  that  if  the  Beurling-Landau 
upper  density  satisfies  D+( A)  <  1,  then  A  is  a  set  of  interpolation. 

The  canonical  case  is  when  LI  =  2n  and  A  =  Z.  Since  {eint}  in  an  o.n.  basis  for  L2[—ir,ir],  it  follows 
from  Parseval  that  A  is  both  a  set  of  sampling  and  a  set  of  interpolation.  This  scales  by  a  change  of 
variable,  and  so  A  =  ^Z  is  both  a  set  of  sampling  and  a  set  of  interpolation  for  PW^ttQ.  Moreover,  general 
lattices  can  be  compared  to  the  canonical  results  as  follows.  If  A  is  a  set  of  sampling  for  PW2„-q,  then 
A  is  everywhere  at  least  as  dense  as  the  lattice  jjZ.  If  A  is  a  set  of  interpolation  for  PW^-n-Q,  then  A  is 
everywhere  at  least  as  sparse  as  the  lattice  ^Z. 

This  generalizes  to  Rd.  Let  D-p  be  a  hyper-rectangle  with  side  lengths  Ll.  If  we  normalize  the  density 
of  Zrf  to  be  one,  then  the  density  of  the  canonical  lattice  for  PW^np  is  l/(27r)rf  times  the  volume  of  the 
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spectrum  fl-p.  Then,  if  A  is  a  set  of  sampling  for  then  A  is  everywhere  at  least  as  dense  as  the 

lattice  If  A  is  a  set  of  interpolation  for  F>W27rnp>  then  A  is  everywhere  at  least  as  sparse  as  the 

lattice 


6.1  Voronoi  Cells  for  Euclidean  Space 

We  use  our  sampling  lattices  to  develop  Voronoi  cells  corresponding  to  the  sampling  lattice.  These  cells 
will  be,  in  the  Euclidean  case,  our  Nyquist  tiles. 

Definition  7  (Voronoi  Cells  in  Rd)  Let  A  =  {A&  G  Rd  :  k  G  N}  be  a  discrete  set  in  Rd.  Then,  the 
Voronoi  cells  {‘hfc},  the  Voronoi  partition  VP( A),  and  partition  norm  ||VT’(A)||  corresponding  to  this  set 
are  defined  as  follows.  Here,  dist  is  the  Euclidean  distance. 

1. )  The  Voronoi  cells  =  {tu  G  :  dist(w,  A*,)  <  inf^  dist(cu,  Aj)}, 

2. )  The  Voronoi  partition  VV( A)  =  {<hfc  G  Rd}keZd, 

3. )  The  partition  norm  | j  ViP (A)  1 1  =  sup^g^dsup^  l/e$k  dist(w,  v). 

Given  /,  /  G  L'2(Rd )  such  that  /  G  PWfjp,  if  the  signal  is  sampled  on  a  lattice  exactly  at  Nyquist, 
we  get  a  sampling  grid  A  =  {A^  G  Rd}k£Zd  that  is  both  a  sampling  set  and  a  set  of  interpolation.  The 
Beurling-Landau  lower  density  and  the  Beurling-Landau  upper  density  are  equal  for  A.  The  dual  lattice 
A*1  in  frequency  space  can  be  used  to  create  Voronoi  cells  {<hfc},  a  Voronoi  partition  W(A-L),  and  partition 
norm  ||V'P(A-L)||.  If  we  sample  on  a  lattice  exactly  at  Nyquist,  each  sample  point  will  correspond  to  an 
element  in  the  dual  lattice  which  is  at  the  center  of  a  Nyquist  tile  NT(/)  for  /.  The  set  of  Nyquist  tiles 
will  cover  Rd.  If,  however,  we  develop  the  Voronoi  cells  {$*,}  for  A^,  we  get  VV(A±)  =  {<&k  G  Rd}k£Zd 
such  that  for  all  k ,  =  {w  G  :  dist(u;,A^)  <  mij^k  dist(cu,  Aj-)}.  But  this  puts  Xk  in  the  center  of 
the  cell.  Then,  if  we  construct  the  Voronoi  cell  containing  this  point,  we  will  get,  up  to  the  boundary, 
the  exact  Nyquist  tile  corresponding  to  this  point.  Nyquist  tiles  are  “half  open,  half  closed.”  If  we  shift 
a  Nyquist  tile  so  that  one  vertex  is  at  the  origin,  we  include  all  of  the  boundaries  that  contain  the  origin, 
and  exclude  the  other  boundaries.  To  get  the  exact  correspondence  between  NT(/)fc  and  &k,  we  make 
Tfc  “half  open,  half  closed,”  and  denote  it  as  $&.  We  denote  the  adjusted  Voronoi  partition  as  VP. 

Theorem  8  (Nyquist  Tiling  for  Euclidean  Space)  Let  f  be  a  non-trivial  function  in  PWuP,  and 
let  A  =  {Afc  G  Rd}k<zz<i  be  the  sampling  grid  which  samples  f  exactly  at  Nyquist.  Let  A1-  be  the  dual  lattice 
in  frequency  space.  Then  the  adjusted  Voronoi  partition  VT’(AJ-)  =  G  Rd}k&zd  equals  the  Nyquist 
Tiling,  i.e., 

{Tfc  G  Rd}keZd  =  {NT (f)k}keZd  . 

Moreover,  the  partition  norm  equals  the  volume  of  A2-,  i.e., 

||V'P(AJ-)||  =  supfceZdSup^e^  dist(w,  v)  =  vol(A±) , 

and  the  sampling  group  G  is  exactly  the  group  of  motions  that  preserve  A1- . 

This  connects,  in  the  Euclidean  case,  sampling  theory  with  the  geometry  of  the  signal  and  its  domain. 
Given  a  function  /  G  PWq,  sampling  of  such  a  function  is  the  process  of  tiling  the  frequency  domain 
by  translated  identical  copies  of  the  parallelepiped  of  minimal  area,  the  Nyquist  Tile,  which  contains 
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the  frequency  support  of  /.  The  relation  between  the  geometry  and  sampling  problem  in  the  Euclidean 
case  is  as  follows:  the  set  of  the  corresponding  translations  -  the  Sampling  Group  -  forms  a  symmetry 
group.  The  corresponding  sampling  set,  which  is  simply  the  annihilator  of  the  sampling  group,  is  also  a 
symmetry  group  of  translations  on  Mrf.  The  set  of  copies  of  the  Nyquist  tile,  obtained  by  applying  the 
sampling  group,  is  the  Nyquist  Tiling. 

The  situation  is  considerably  different  when  the  underlying  space  is  not  Euclidean.  We  quickly  get 
into  open  problems.  Theorem  8  gives  an  approach  for  solving  the  problem  in  non-Euclidean  spaces.  We 
suggest  using  the  two  tools  we  just  established  -  the  Beurling-Landau  density  and  Voronoi  cells.  Our 
next  section  discusses  the  geometry  of  orientable  surfaces.  In  particular,  it  provides  insight  into  why  a 
focus  on  Euclidean,  spherical,  and  especially,  hyperbolic  geometries  is  important. 

6.2  Sampling  in  Hyperbolic  Space 

We  begin  by  stating  the  Fourier  transform,  its  inversion,  and  the  Plancherel  formula  for  hyperbolic  space 
[HelOO]. 

Let  dz  denote  the  area  measure  on  the  unit  disc  D  =  {z  \  \z\  <  1,  and  let  the  measure  dv  be  given  by 
then  the  SU(1,  l)-invariant  measure  on  D  is  given  by  dv(z )  =  dz/(  1  —  |z|2)2.  For  functions  /  £  L1( D,  dv) 
the  Helg  as  on- Fourier  transform  is  defined  as 

/(A,6)=  [  f(z)e^~iX+1^  dv(z) 

J  D 

for  A  >  0  and  b  £  T.  Here  (z,  b)  denotes  the  minimal  hyperbolic  distance  from  the  origin  to  the 
horocycle  through  z  and  a  point  b  £  <9D.  The  mapping  /  i— >  /  extends  to  an  isometry  L2(D,du)  — ► 
L2(M+  x  T,  (27r)_1Atanh(A7r/2)dAd6),  i.e.,  the  Plancherel  formula  becomes 

[  \f(z)\2n  d?\2\2  =  S~  f  |/(A,6)|2Atanh(Avr/2)dAd6. 

J D  (1  -\A  )  2n  J R+xT 

Here  db  denotes  the  normalized  measure  on  the  circle  T,  such  that  fT  db  =  1,  and  dX  is  Lebesgue  measure 
on  M.  The  Helg  as  on- Fourier  inversion  formula  is 

f(z)  =  —  f  [  /(A,  6)e(a+1^z,^Atanh(A7r/2)  dX  db . 

27r  Jr+  J t 

A  function  /  £  L2(0,dv)  is  called  band-limited  if  its  Helgason-Fourier  transform  /  is  supported  inside  a 
bounded  subset  [0,  fl]  of  M+.  The  collection  of  band-limited  functions  with  band-limit  inside  a  set  [0,  H] 
will  be  denoted  PWq  =  PWn(B).  This  definition  of  band-limit  coincides  with  the  definitions  given  in 
[FP04]  and  [( '013]  which  both  show  that  sampling  is  possible  for  band-limited  functions.  The  Laplacian 
on  D  is  symmetric  and  given  by 

and  we  note  that 

A/(A,6)  =  -(A2  +  1)/(A,6). 

Therefore,  if  /  £  PWq(D),  we  see  that  the  following  Bernstein  inequality  is  satisfied  ||An/||  <  (1  + 

|si|2)”/2ll/l!  ' 

In  the  following  section  we  will  describe  sampling  results  for  band-limited  functions  on  hyperbolic 
space,  which,  it  must  be  stressed,  do  not  deal  with  optimal  densities. 
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6.2.1  Sampling  via  Operator  Theory  in  B 

The  work  in  [FP04]  defines  band-limits  using  the  spectrum  of  the  Laplacian  on  a  manifold,  while  [ 
builds  on  representation  theory  which  for  the  case  at  hand  gives  the  explicit  form  of  the  Fourier  transform 
on  B  as  defined  above.  We  also  refer  to  the  paper  [FP05]  which  provides  the  same  results  in  the  setting 
of  the  upper  half  plane  (which  is  bi-holomorphically  equivalent  to  B).  These  papers  build  on  Neumann 
series  for  an  operator  based  on  sampling  as  well  as  the  Bernstein  inequality.  The  sampling  operators  have 
previously  been  explored  in  [Gro91,  Gro92]. 

According  to  Pesenson  [PesOO]  there  is  a  natural  number  N  such  that  for  any  sufficiently  small  r  there 
are  points  xj  S  B  for  which  B(xj,r/ 4)  are  disjoint,  B(xj,r/ 2)  cover  B  and  1  <  JG  XB(xj,r)  —  N-  Such  a 
collection  of  {xj}  will  be  called  an  (r,  iV)-lattice. 

Let  cf)j  be  smooth  non-negative  functions  which  are  supported  in  B[xj,r/ 2)  and  satisfy  that  JG  <t>j  = 

Id  and  define  the  operator  Tf(x)  =  Pq  (Ylj  f(xj)^j(x)^  >  where  Pq  is  the  orthogonal  projection  from 

L2(D,du)  onto  PWq(B).  By  decreasing  r  (and  thus  choosing  Xj  closer)  one  can  obtain  the  inequality 
|| I  —  T\\  <  1,  in  which  case  T  can  be  inverted  by 


OO 

T~lf  =  YJ{I~T)kf. 

k= 0 

For  given  samples  we  can  calculate  Tf  and  the  Neumann  series  provides  the  recursion  formula  fn+i  = 
fn  +  T f  —  T fn  and  then  lim^^oo  fn  =  f  with  norm  convergence.  The  rate  of  convergence  is  determined 
by  the  estimate  ||/n  -  /||  <  ||J  -  T||n+1||/||. 

The  paper  [FP05]  further  provides  a  necessary  condition  for  the  set  {xi\  to  be  a  sampling  set.  They 
find  that  there  is  a  constant  C  which  is  determined  by  the  geometry  of  B,  such  that  if  r  <  C^1(l  + 
|j^|2)fc/2)-1  for  any  then  any  (IV,  r)-lattice  {xi}  is  a  sampling  set.  The  paper  [C013]  obtains 

similar  results,  but  removes  some  restrictions  on  the  functions  (f>j.  In  particular  the  partitions  of  unity 
do  not  need  to  be  smooth  and  can  actually  be  chosen  as  characteristic  functions  (f>j  =  xUj  f°r  a  cover  of 
disjoint  sets  Uj  contained  in  the  balls  B{xj,r/ 2).  This  is  done  by  lifting  the  functions  to  the  group  of 
isometries  (which  in  this  case  is  SXJ ( 1 , 1 ) ) ,  and  by  estimating  local  oscillations  using  Sobolev  norms  for 
left-invariant  vector  fields  on  this  group. 

6.2.2  Beurling  density  for  Bergman  spaces 

In  this  section  we  describe  a  collection  of  celebrated  sampling  theorems  for  Bergman  spaces  on  the  unit 
disc  by  [Sei93]  and  [Sch97].  Let  77(B)  be  the  space  of  holomorphic  functions  on  B.  Let  1  <  p  <  oo 
be  given,  and  equip  the  unit  disc  B  with  normalized  area  measure  dcr(z).  We  define  the  Bergman 
space  AP(B)  =  Lp(B,du)  n  77(B).  This  is  a  reproducing  kernel  Banach  space  with  reproducing  kernel 
K(z,w )  =  =rp.  By  [Sci93]  and  [Sch97]  sampling  and  interpolation  sets  for  AP(B)  are  characterized  by 
the  upper  and  lower  Beurling  densities 

D+(Z )  =  limsup  sup  D((/)W(Z),  r) ,  D~(Z)  =  liminf  inf  D(cf>w(Z),r) . 

1 — >1  tcSD  1  >1  W^D 

Here  <j>w(z)  =  and  D(Z,r)  =  (Z)|,fc|<r  log  |zfc|)/(log(l  -r)).  Let  p(z,w)  =  \<j>w(z)\  be  the  pseudo- 
hyperbolic  distance  from  z  to  w,  then  a  sequence  Z  =  { Zi }  is  called  uniformly  discrete  if  there  is  a  6  >  0 
such  that  p(zi,  Zj)  >5  for  i  ^  j. 

Theorem  9  Let  A.  be  a  set  of  distinct  points  in  B. 
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1. )  A  sequence  A  is  a  set  of  sampling  for  Ap  if  and  only  if  it  is  a  finite  union  of  uniformly  discrete  sets 

and  it  contains  a  uniformly  discrete  subsequence  A'  for  which  D~(A')  >  1/p. 

2. )  A  sequence  A  is  a  set  of  interpolation  for  Ap  if  and  only  if  it  is  uniformly  discrete  and  D+( A)  <  1/p. 

This  results  show  there  can  be  no  Nyquist  density  for  the  Bergman  spaces,  since  the  sampling  sets 
are  always  sharply  separated  from  the  interpolating  sets.  We  note  that  the  results  of  Seip  and  Schuster 
are  for  a  particular  class  of  holomorphic  functions,  to  which  the  band-limited  functions  PW^m)  do  not 
belong.  It  is  an  open  question  whether  it  is  possible  to  establish  a  Nyquist  density  for  band-limited 
functions  on  D  and  to  use  this  information  to  create  regular  lattices  and  dual  lattices  determined  by  the 
size  of  the  band-limit  Q. 


6.2.3  Voronoi  Cells  and  Beurling-Landau  Density  for  D 


We  develop  our  model  for  hyperbolic  space  on  the  Poincare  disk  ID).  The  motions  that  preserve  lengths 
in  hyperbolic  geometry  are  Mobius-Blaschke  maps.  Geodesics  are  subarcs  of  paths  that  intersect  dD  at 
right  angles.  Let  T  be  a  smooth  path  in  the  unit  disk  D.  The  hyperbolic  length  of  T  is  £#(T)  =  jf,  • 

Let  a£D,  and  let  =  e™  (a  Mobius-Blaschke  transformation  of  D  onto  D).  Then  p>e,a  preserves 
the  hyperbolic  length,  i.e,  =  £#(T)  .  The  hyperbolic  distance  p  between  two  points  z\,  Z2 

in  D  is 


P(zi,z2) 


2  arctanh 


(  1 31  ~  Z2l  \ 


log 


\zi~z2\  \ 
|1~ Z2Zl|  I 

I Zl  Z2 |  I 

|l-Z2Zl|  / 


The  distance  p  will  be  used  to  determine  distance  for  the  sampling  lattice  A.  Note  that,  because  we 
cannot  establish  the  Beurling-Landau  densities,  we  can  not  create  regular  lattices  and  dual  lattices. 

The  Helgason-Fourier  transform  maps  L2(D)  to 


L2(R+  x  T, 


— Atanh(A7r/2)dA  db) , 

Z7T 


which  is  isomorphic  to  the  space  of  L2(T)-vector  valued  square  integrable  functions  with  measure  Atanh(A7r/2)dA, 
in  short  denoted  by 


L2(M+;L2(T),Atanh(A7r/2)dA)  . 


The  negative  Laplacian  —A  is  positive  with  spectrum  M+,  and  therefore  we  define  Voronoi  cells  based  on 
a  distance  on  M+.  This  distance  is  denoted  dist,  and  it  is  an  open  question  in  which  manner  it  is  related 
to  the  measure  Atanh(A7r/2)  dX.  With  an  appropriate  distance  function  dist,  we  can  define  the  following. 


Definition  8  (Voronoi  Cells  in  D)  Let  A  =  {A^,  £  D  =  M+  x  T  :  k  £  N},  be  a  discrete  set  in  fre¬ 
quency  space.  Then,  the  Voronoi  cells  {‘hfc},  the  Voronoi  partition  VP(A),  and  partition  norm  || VW>(A)  || 
corresponding  to  this  set  are  defined  as  follows. 


1. )  The  Voronoi  cells  =  {u>  £  D  :  dist(<u,Afc)  <  inf^j,  dist(cu,  Aj)}; 

2. )  The  Voronoi  partition  VV{ A)  =  {<!>£  C 

3. )  The  partition  norm  | j  ViP (A)  1 1  =  sup^.g^dsup^ dist (cu,  u). 


A  crucial  step  in  answering  the  question  of  Nyquist  density  using  Voronoi  cells  is  to  determine  an 
appropriate  candidate  for  the  distance  on  D. 
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6.3  Sampling  on  the  Sphere 


One  perspective  of  Fourier  analysis  is  to  think  of  it  as  a  systematic  use  of  symmetry  to  simplify  and 
understand  linear  operators.  The  unit  sphere  §2  admits  the  special  orthogonal  group  of  three  variables, 
SO( 3)  -  proper  rotations  of  M3  about  the  origin  -  as  a  transitive  group  of  symmetries.  Fourier  analysis 
on  §2  amounts  to  the  decomposition  of  L2(§2)  into  minimal  subspaces  invariant  under  all  rotations  in 
SO(3).  The  rotations  of  the  sphere  induce  operators  on  functions  by  rotating  the  graphs  over  §2.  The 
Hilbert  space  L2(§2)  is  defined  with  the  usual  inner  product,  using  the  rotation-invariant  area  element 
p.  Background  for  this  section  can  be  found  in  Driscoll  and  Healy  [  ],  Keiner,  Kunis,  and  Potts 

[KKP07],  and  McEwen  and  Wiaux  [MW11]. 

Band- limited  functions  on  the  sphere  are  spherical  polynomials.  The  corresponding  sampling  problem 
is  the  computation  of  Fourier  coefficients  of  a  function  from  sampled  values  at  scattered  nodes.  If  we 
consider  the  problem  of  reconstructing  a  spherical  polynomial  of  degree  JVgN  from  sample  values,  one 
might  set  up  a  linear  system  of  equations  with  M  =  ( N  +  l)2  interpolation  constraints  which  has  to 
be  solved  for  the  unknown  vector  of  Fourier  coefficients  f  G  C(Ar+1)  .  If  the  nodes  for  interpolation  are 
chosen  such  that  the  interpolation  problem  always  has  a  unique  solution,  the  sampling  set  is  called  a 
fundamental  system. 

Let  §2  =  {x  Gl3  :  ||x||2  =  1}  be  the  two-dimensional  unit  sphere  embedded  in  M3.  A  point  p  G  S2  is 
identified  in  spherical  coordinates  by  rj  =  (sin(0)  cos(0),  sin($)  sin(iji),  cos(0))T,  where  the  angles  ( 9 ,  <f> )  are 
the  co-latitude  and  longitude  of  g.  Topologically,  S2  =  C.  Geodesics  are  great  circles,  and  the  geodesic 
distance  can  be  most  directly  written  as  dist(p,^)  =  arccos (r/  •  £)  .  For  r/,£  G  §2,  \\g  —  £||2  =  2  —  2 (77  •  £). 
The  distance  to  the  “north  pole”  n  =  (0,  0, 1)T  of  §2  is  arccos (77  •  to)  =  9. 

The  spherical  harmonics  YfJ  form  an  o.n.  basis  for  L2(§2)  We  can  define  them  as  follows.  The  Legendre 
polynomials  P \  :  [—1, 1]  — >  M  are  generated  by  applying  the  Gram-Schmidt  method  to  {xfc}£T0.  They 
are  given  by  the  Rodrigues  formula  Pk{t)  =  1/(2 kk\)dk / dtk{t2  —  l)fc.  The  associated  Legendre  functions 
are  defined  by 


Pkit)  = 


(TTT)!<f  _1)  d^Pt(t)- 


The  spherical  harmonics  Yff  :  §2  — »  C  of  degree  k  G  N  U  {0}  and  order  n  G  Z,  |n|  <  k,  are  the  functions 
rfen(p)  =  Y£(6,  (/>)  =  .  We  have  that 


/*27T  /*7T 

/  /  Y£(9,  <f>)Yr(9,  <f)  sin (9)  d9  def  =  , 

Jo  Jo 

i.e. ,  Yu  form  an  o.n.  basis  for  L2(§2).  We  say  that  /  is  a  spherical  polynomial  of  degree  N  if  f(9,<j) )  = 
J2k= 0  Yln=-k  fk^k-  r^'^ie  sPace  °f  spherical  polynomials  of  degree  at  most  N  has  dimension  ( N  +  l)2. 

The  Fourier  transform  is  the  spherical  Fourier  matrix  f(9,  <f)  =  ^2n=-k  fk^k  >  whh  coefficients 

given  by  fJJ  =  fs2  fYff  dp .  The  dual  space  of  L2(§2)  is  discrete.  The  inverse  Fourier  transform  is  the 
construction  of  a  spherical  polynomial  from  the  coefficients.  The  function  f  is  N  band-limited  (IV  G  N) 
if  fk  =  0  f°r  k  >  N.  Thus,  f(9,  (f)  =  Ylk= 0  Yln=-k  fk^k-  -^or  problem  of  solving  for  a  spherical 
polynomial  f  of  degree  N  from  sample  values,  we  are  looking  to  solve  for  the  unknown  Fourier  coefficients 

{fjf}  = 

Let  A  =  {A k\i=i  be  a  sampling  set  on  §2.  The  mesh  norm  and  the  separation  distance  q,\  are 
defined  by  5\  =  2  m ax^ e g 2 m i n ^  =i , m  dist (77 ,  A k)  ,  q\  =  minj^.  dist(Aj,  A k)  ■  A  sampling  set  A  is  called  5 
dense  if  for  some  0  <  5  <  2i r,  5\  <  5,  and  called  q  separated  if  there  exists  0  <  q  <  2ir  such  that  q\  >  q. 
We  assume  that  our  sampling  set  is  separated.  Finally,  a  sampling  set  is  called  quasi-uniform  if  there 
exist  a  constant  C  independent  of  the  number  on  sample  points  M  such  that  5\  <  C  q\. 
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Sampling  on  the  sphere  is  how  to  sample  a  band-limited  function,  an  IVth  degree  spherical  polynomial, 
at  a  finite  number  of  locations,  such  that  all  of  the  information  content  of  the  continuous  function  is 
captured.  Since  the  frequency  domain  of  a  function  on  the  sphere  is  discrete,  the  spherical  harmonic 
coefficients  describe  the  continuous  function  exactly.  A  sampling  theorem  thus  describes  how  to  exactly 
recover  the  spherical  harmonic  coefficients  of  the  continuous  function  from  its  samples.  Given  A,  the 
spherical  Fourier  transform  matrix  is  Y  =  ...:M;k=o,...,N;\n\<k  •  Let  Yh  denote  its  complex 

conjugate  transpose.  The  inverse  Fourier  transform  matrix  is  the  construction  of  a  spherical  polynomial 
of  degree  N  from  given  data  points  (A j,yj)  £  §2  x  C  such  that  the  identity  f(Xj)  =  y3  is  solved.  This 
is  solving  the  linear  system  Yf  =  y,  y  =  (yi,y2,  ■  ■  ■  ,Vm)  for  the  vector  of  Fourier  coefficients  f  =  {/£■} 
of  the  spherical  polynomial.  Essentially,  it  is  the  inverse  problem  to  /  =  Yf,  which  corresponds  to 
evaluating  a  spherical  polynomial  on  A. 

The  open  question  again  is  the  establishment  of  the  optimal  Beurling-Landau  densities.  This  leads  to 
questions  about  sphere  tiling.  The  papers  [KKP07]  and  [  ]  address  the  problem  of  finding  optimal 

sampling  lattices. 

6.4  General  Analytic  Surfaces 

A  surface  is  a  generalization  of  Euclidean  space.  From  the  viewpoint  of  harmonic  analysis,  there  is  a 
natural  interest  in  both  the  theory  and  applications  of  the  study  of  integrable  and  square  integrable 
functions  on  surfaces.  Background  material  for  this  subsection  can  be  found  in  Ahlfors  [Ahl78,  AhllO], 
Farkas  and  Kra  [FK80],  Forster  [For81],  Lee  [Lec97],  and  Singer  and  Thorpe  [ST67] . 

We  assume  our  surfaces  are  connected  and  orientable.  Therefore,  we  can  choose  a  coordinate  system 
so  that  differential  forms  are  positive  [  ].  We  consider  Riemann  surfaces,  but  our  discussion  carries 

through  to  connected  and  orientable  Riemannian  manifolds  of  dimension  two  [  ].  Riemann  surfaces 

allow  us  to  discuss  the  Uniformization  Theorem,  which  gives  that  all  orientable  surfaces  inherit  their 
intrinsic  geometry  from  their  universal  coverings.  There  are  only  three  universal  covers  -  the  plane  C 
(Euclidean  geometry) ,  the  Riemann  sphere  C  (spherical  geometry) ,  and  the  hyperbolic  disk  D  (hyperbolic 
geometry) . 

Given  connected  Riemann  surface  S  and  its  universal  covering  space  S,  S  is  isomorphic  to  S/G, 
where  the  group  G  is  isomorphic  to  the  fundamental  group  of  S,  7Ti(«S)  (see  [For81],  Section  5).  The 
corresponding  universal  covering  is  simply  the  quotient  map  which  sends  every  point  of  S  to  its  orbit 
under  G.  Thus,  the  fundamental  group  of  S  determines  its  universal  cover.  Moreover,  the  universal 
covering  is  indeed  the  “biggest”  smooth  unlimited  covering  of  a  connected  Riemann  surface,  in  the  sense 
that  all  other  unramified  unlimited  covering  space  of  a  Riemann  surface  can  be  covered  unlimitedly  and 
without  ramification  by  the  universal  covering  of  this  surface. 

The  Uniformization  Theorem  allows  us  to  classify  all  universal  covers  of  all  Riemann  surfaces.  This 
in  turn  allows  us  to  understand  the  geometry  of  every  Riemann  surface. 

Theorem  10  (The  Uniformization  Theorem)  Let  S  be  a  Riemann  surface. 

1. )  Every  surface  admits  a  Riemannian  metric  of  constant  Gaussian  curvature  k. 

2. )  Every  simply  connected  Riemann  surface  is  conformally  equivalent  to  one  of  the  following: 

a. )  C  -  Euclidean  Geometry  -  k  =  0  -  with  isometries  +  a |  ,  0  e  [0,  27t)  , 

b. )  C  -  Spherical  Geometry  -  k  =  1  -  with  isometries  l  ,  ,  |a|2  +  |/3|2  =  1 , 
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c.) 


O  -  Hyperbolic  Geometry  -  k  =  —  1  -  with  isometries 


Given  connected  Riemann  surface  S  and  its  universal  covering  space  S ,  S  is  isomorphic  to  S/G, 
where  the  group  G  is  isomorphic  to  the  fundamental  group  of  S,  717  (5)  (see  [ForSl],  Section  5).  The 
corresponding  universal  covering  is  simply  the  quotient  map  which  sends  every  point  of  S  to  its  orbit  under 
G.  Forster  [ForSl]  (Section  27)  gives  the  consequences  of  the  Uniformization  Theorem  very  succinctly. 
The  only  covering  surface  of  Riemann  sphere  C  is  itself,  with  the  covering  map  being  the  identity.  The 
plane  C  is  the  universal  covering  space  of  itself,  the  once  punctured  plane  C  \  {zo}  (with  covering  map 
exp(z  —  zo));  and  all  tori  C/F,  where  T  is  a  parallelogram  generated  by  z  i — >  z  +  rv yi  +  77172  ,  n,  m  £  Z 
and  71,72  are  two  fixed  complex  numbers  linearly  independent  over  M.  The  universal  covering  space  of 
every  other  Riemann  surface  is  the  hyperbolic  disk  O.  Therefore,  the  establishment  of  exact  the  Beurling- 
Landau  densities  for  functions  in  Paley-Wiener  spaces  in  spherical  and  especially  hyperbolic  geometries 
will  allow  the  development  of  sampling  schemes  on  arbitrary  Riemann  surfaces. 


6.5  The  Kunze-Stein  phenomenon 

In  this  subsection  we  will  point  out  that  for  some  amenable  subgroups  of  semi-simple  Lie  groups,  the 
continuity  of  convolution  operators  follows  from  the  famous  Kunze-Stein  phenomenon  [KS60,  Cow78, 
Cow08].  This  work  is  joint  with  Jens  Christensen  [(  '('015]. 

Let  G  be  a  semi-simple,  connected  and  non-compact  group  with  Haar  measure  dx ,  then 

Theorem  11  (Cowling)  Suppose  that  1  <  r  <  00,  1/r-f  1/r'  =  1  and  f  £  Lr(G ).  Then  the  mapping 
9  ^  9  *  f  bounded  from  LP(G )  to  Lq(G)  provided  that  one  of  the  following  conditions  hold: 

1.  if  r  =  1  and  1  <  p  =  q  <  00 

2.  if  1  <  r  <  2 ,q  >r,p<r',  0  <  1/p  —  1/q  <  1/r'  and  ( p ,  q)  ^  (r,  r)  and  ( p ,  q )  ^  (r',  r'). 

3.  if  2  <  r  <  00,  q  >  r,  p  <  r' ,  0  <  1/p  —  1/q  <  1/r'  and  ( p ,  q)  ^  (r,  r') 

4 ■  if  r  =  00, p  =  1  and  q  =  00. 


We  will  only  focus  on  the  second  condition  with  p  =  q,  in  which  case  the  condition  simplifies  to:  Lp*Lr  C 
IT  if  1  <  r  <  2  and  r  <  p  <  r/(r  —  1). 

Since  G  is  semi-simple  there  is  a  compact  subgroup  I\ ,  an  abelian  subgroup  A  and  a  nilpotent  subgroup 
N  such  that  A  x  N  x  K  3  (a,  n,  k )  1— >  ank  E  G  is  a  diffeomorphism  and  with  proper  normalizations  of 
Haar  measures 


/  f(x)dx=  /  /  /  f  (ank)  da  dn  dk. 

Jg  Jk  J  a  Jn 


If  /  is  a  Jv-right  invariant  function  on  G,  then  it  can  be  identified  with  a  function  /  on  S  =  AN  and 


IG 


f{x)dx=  I  f(s)ds=  I  f  f(an)dadn. 

J  s  J  A  Jn 


If  /  is  A"-bi-invariant  and  g  is  /\-right-invariant,  then  the  convolution  of  the  functions  g  and  /  on  S  =  AN 
can  be  written  as  a  convolution  of  the  corresponding  functions  on  G:  g  *  f  =  g  *  f . 

From  the  Kunze-Stein  phenomenon  we  can  now  obtain  continuity  of  convolution  operators  on  the 
amenable  (subgroup)  S  =  AN.  This  is  the  second  contribution  of  this  paper: 
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Theorem  12  Assume  that  f  is  a  function  on  S  =  AN  and  that  there  is  a  K -bi-invariant  function  h  on 
G  for  which  |/(s)|  <  h(s)  for  all  s  E  5.  If  h  is  in  Lr(G)  for  some  r  e  (1,  2)  then  g  i— >  g  *  /  is  continuous 
on  LP(S)  if  r  <  p  <  r /{r  —  1). 

Example  1  This  example  shows  that  when  G  =  SU(n,  1)  convolution  by  a  certain  non-lfiljG)  function  is 
continuous  on  LP(G).  This  is  related  to  the  Bergman  projection,  and  we  will  follow  this  path  to  settling 
open  questions  surrounding  atomic  decompositions  and  the  Bergman  projection  on  unbounded  symmetric 
domains  [GR80,  BBGROf]. 

Any  x  £  G  can  be  written  as  a  matrix 


where  a  is  an  n  x  n-matrix,  b,  c  E  Cn  and  d  £  C. 

The  (multivalued)  function 

f(x)  =  1/eT, 

corresponds  to  a  single-valued  function  f  when  restricted  to  the  simply  connected  subgroup  S.  If  we 
define  h(x)  =  \f(x)\,  then  h  is  K -bi-invariant  and  |/(s)|  <  h(s).  Also,  the  function  h  is  in  Lr(G)  for 
r  >  2n/a,  and  therefore  convolution  by  f  is  continuous  on  LP(S )  when  2n/a  <  p  <  2n/(2n  —  a)  and 
n  <  a  <  2 n.  Note,  that  if  n  =  1  this  is  exactly  the  statement  from  Example  ??7  since  SL(2,  M)  and 
SU(1, 1)  are  isomorphic.  Also,  if  a  >  2 n,  then  h  is  in  T1(G)  and  the  continuity  of  the  convolution 
operator  is  immediate.  One  interesting  aspect  of  deriving  these  convolution  continuities  on  SU(n,  1)  is 
their  relation  to  the  Bergman  projection  on  the  unit  ball  in  Cn.  In  particular  they  can  be  used  (when 
applied  to  weighted  Lp(S)-spaces)  to  prove  the  following  classical  theorem  [Zhu05,  Thm.  2.10] 

Theorem  13  Fix  two  real  parameters  a  and  b  and  define  the  integral  operator  S  by 

S  f(z)  =  (1  —  \z\2)a  [  — — ~ r - 1  I  ,  f(w)  dv(w), 

^  JMn  |1  -  (Z,w)\n+1  +  a+b  ^  h 

where  dv  is  the  volume  measure  on  the  unit  ball  Bn  in  Cn  and  ( z,w )  =  z\W\  +  •  •  •  +  znwn.  Then,  for  t 
real  and  1  <  p  <  oo,  S  is  bounded  on  Lf(Bn)  if  and  only  if  —pa  <  t  +  1  <  p(b  +  1). 

We  are  currently  working  out  the  details  of  this  argument  and  will  publish  it  in  a  forthcoming  paper. 


7  Analysis  of  Point  Processes 

The  analysis  of  periodic  processes  is  an  important  area  of  signal  analysis.  A  probabilistic  view  of  the 
Riemann  zeta  function,  algorithms  on  generalizations  of  Euclidean  domains,  and  variations  on  equidis- 
tribution  theory  have  led  to  algorithms  for  several  classes  of  problems  in  periodic  parameter  estimation. 
These  algorithms  are  general  and  very  efficient,  and  can  be  applied  to  the  analysis  of  periodic  processes 
with  a  single  generator  and  the  deinterleaving  and  analysis  of  processes  with  multiple  generators. 

We  divide  our  analysis  into  two  cases  periodic  processes  created  by  a  single  source,  and  those  processes 
created  by  several  sources.  We  wish  to  extract  the  fundamental  period  of  the  generators,  and,  in  the 
second  case,  to  deinterleave  the  processes.  We  first  present  very  efficient  algorithm  for  extracting  the 
fundamental  period  from  a  set  of  sparse  and  noisy  observations  of  a  single  source  periodic  process.  The 
procedure  is  computationally  straightforward,  stable  with  respect  to  noise  and  converges  quickly.  Its  use 
is  justified  by  a  theorem  which  shows  that  for  a  set  of  randomly  chosen  positive  integers,  the  probability 
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that  they  do  not  all  share  a  common  prime  factor  approaches  one  quickly  as  the  cardinality  of  the  set 
increases.  The  proof  of  this  theorem  rests  on  a  probabilistic  interpretation  of  the  Riemann  zeta  function. 
We  then  build  upon  this  procedure  to  deinterleave  and  then  analyze  data  from  multiple  periodic  processes. 
This  relies  both  on  the  probabilistic  interpretation  of  the  Riemann  zeta  function,  the  equidistribution 
theorem  of  Weyl,  and  Wieners  periodogram. 

7.1  Single  Generator:  The  ME  A  Algorithm 

We  have  developed  a  collection  of  algorithms  that  analyze  periodic  phenomena  generated  by  a  single 
generator.  These  work  even  when  the  data  is  extremely  sparse  and  noisy.  The  algorithms  use  number 
theory  in  novel  ways  to  extract  the  underlying  period  by  modifying  the  Euclidean  algorithm  to  determine 
the  period  from  a  sparse  set  of  noisy  measurements  [CS13,  SCOO,  SC98,  SC96,  Cas97,  CS96,  CS96].  The 
elements  of  the  set  are  the  noisy  occurrence  times  of  a  periodic  event  with  (perhaps  very  many)  missing 
measurements.  The  proposed  algorithms  are  computationally  straightforward  and  converge  quickly.  A 
robust  version  is  developed  that  is  stable  despite  the  presence  of  arbitrary  outliers.  The  Euclidean 
algorithm  approach  is  justified  by  a  theorem  which  shows  that,  for  a  set  of  randomly  chosen  positive 
integers,  the  probability  that  they  do  not  all  share  a  common  prime  factor  approaches  one  quickly  as  the 
cardinality  of  the  set  increases.  The  theorem  is  in  essence  a  probabilistic  interpretation  of  the  Riemann 
Zeta  Function.  In  the  noise-free  case  this  implies  convergence  with  only  ten  data  samples,  independent  of 
the  percentage  of  missing  measurements.  In  the  case  of  noisy  data  simulation  results  show,  for  example, 
good  estimation  of  the  period  from  one  hundred  data  samples  with  fifty  percent  of  the  measurements 
missing  and  twenty  five  percent  of  the  data  samples  being  arbitrary  outliers. 

We  then  use  these  algorithms  in  the  analysis  of  periodic  pulse  trains,  getting  an  estimate  of  the 
underlying  period  [C'S13,  SCOO,  SC98,  SC96,  Cas97,  CS96,  CS96].  This  estimate,  while  not  maximum 
likelihood,  is  used  as  initialization  in  a  three-step  algorithm  that  achieves  the  Cramer-Rao  bound  for 
moderate  noise  levels,  as  shown  by  comparing  Monte  Carlo  results  with  the  Cramer-Rao  bounds.  An 
approach  using  multiple  independent  data  records  is  also  developed  that  overcomes  high  levels  of  con¬ 
tamination.  The  data  sets  arises  in  radar  pulse  repetition  interval  (PRI)  analysis,  in  bit  synchronization 
in  communications,  in  biomedical  applications,  and  other  scenarios.  We  assume  our  data  is  a  finite  set  of 
real  numbers  S  =  {sj}™=l  ,  with  Sj  =  kjr  +  4>  +  rjj,  where  r  (the  period)  is  a  fixed  positive  real  number, 
the  kj's  are  non-repeating  positive  integers,  (j)  (the  phase)  is  a  real  random  variable  uniformly  distributed 
over  the  interval  [0,r),  and  the  ijj' s  are  zero-mean  independent  identically  distributed  (iid)  error  terms. 
We  assume  that  the  rjj's  have  a  symmetric  probability  density  function  (pdf),  and  that  \r)j\  <  |  for  all 
j .  We  develop  an  algorithm  for  isolating  the  period  of  the  process  from  this  set,  which  we  shall  assume 
is  (perhaps  very)  sparse.  In  the  noise-free  case  our  basic  algorithm,  given  below,  is  equivalent  to  the  Eu¬ 
clidean  algorithm  and  converges  with  very  high  probability  given  only  n  =  10  data  samples,  independent 
of  the  number  of  missing  measurements.  We  assume  that  the  original  data  set  is  in  descending  order, 

i.e. ,  Sj  >  Sj. |_i.  Let  r  denote  the  value  the  algorithm  gives  for  r,  and  let  “« — ”  denote  replacement ,  e.g., 
“a  < —  b ”  means  that  the  value  of  the  variable  a  is  to  be  replaced  by  the  current  value  of  the  variable  b. 

The  Modified  Euclidean  Algorithm  (MEA) 

Initialize:  Set  iter  =  0. 

1. )  [Adjoin  0  after  first  iteration.]  If  iter  >  0,  then  S  < —  S  U  {0}. 

2. )  [Form  the  new  set  with  elements  ( Sj  —  Sj+i).]  Set  Sj  < —  ( Sj  —  Sj+i). 

3. )  Sort  the  elements  in  descending  order. 

4. )  [Eliminate  noise.]  If  0  <  Sj  <  rjo,  then  S  < —  S  \  {sj}. 

5. )  The  algorithm  terminates  if  S  has  only  one  element  si.  Declare  r  =  si.  If  not,  then  set 
iter  ^ —  (iter  +  1).  Go  to  (1.). 
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Here,  0  <  r?o  <  r  is  a  noise  threshold.  Noise-free  simulation  examples  demonstrate  successful  estimation 
of  r  for  n  =  10  with  99.99%  of  the  possible  measurements  missing.  In  fact,  with  only  10  data  samples,  it 
is  possible  to  have  the  percentage  of  missing  measurements  arbitrarily  close  to  100%.  There  is,  of  course, 
a  cost,  in  that  the  number  of  iterations  the  algorithm  needs  to  converge  increases  with  the  percentage  of 
missing  measurements.  In  the  presence  of  noise  and  false  data  (outliers),  there  is  a  tradeoff  between  the 
number  of  data  samples,  the  amount  of  noise,  and  the  percentage  of  outliers.  The  algorithm  will  perform 
well  given  low  noise  for  n  =  10,  but  will  degrade  as  noise  is  increased.  However,  given  more  data,  it 
is  possible  to  reduce  noise  effects  and  speed  up  convergence  by  binning  the  data,  and  averaging  across 
bins.  Binning  can  be  effectively  implemented  by  using  an  adaptive  threshhold  with  a  gradient  operator, 
allowing  convergence  in  a  single  iteration  in  many  cases.  Simulation  results  show,  for  example,  good 
estimation  of  the  period  from  one  hundred  data  samples  with  fifty  percent  of  the  measurements  missing 
and  twenty  five  percent  of  the  data  samples  being  arbitrary  outliers  [CS13,  CS96]. 

Our  algorithm  is  based  on  several  theoretical  results,  which  we  now  present.  First,  we  can  modify 
the  basic  Euclidean  algorithm,  allowing  a  reformulation  using  subtraction  rather  than  division,  because 
gcd(rAq, . . . ,  rkn )  =  r  gcd((/t'i  —  kf), . . . ,  (fcn-i  —  kn ),  kn )  .  We  then  show  that  our  procedure  almost  surely 
converges  to  the  period  by  proving  the  following  result.  The  Riemann  Zeta  Function  is  defined  on  the 
complex  half  space  {z  £  C  :  K(z)  >  1}  by  =  Ylri=i  n~z  ■  Euler  demonstrated  the  connection  of 
(  with  number  theory  by  showing  that  C(z)  =  n=i  i-(p  )-z  >  >  1  ■  where  P  =  {pi,P2,P3,  ■  ■  ■}  = 

{2, 3, 5, . . .}  is  the  set  of  all  prime  numbers.  In  the  following,  we  let  P{-}  denote  probability,  card{-}  denote 
the  cardinality  of  the  set  {•},  and  let  {1, . . .  ,£}n  denote  the  sublattice  of  positive  integers  in  Mn  with 
coordinates  c  such  that  1  <  c  <  l.  Therefore,  Nn(£)  =  card{(/ci, . . . ,  kn)  e  {1, . . . ,  £}n  :  gcd(£q, . . . ,  kn)  = 
1}  is  the  number  of  relatively  prime  elements  in  {1, . . .  ,£}n. 

Theorem  14  ([CS13,  CS96])  For  n  >  2,  we  have  that  lim^oo  Nn^  =  [C(n)]  1  •  Therefore,  given 
n  (n  >  2)  randomly  chosen  positive  integers  {ki, . . .  ,kn},  P{gcd(fei, . . . ,  kn)  =  1}  =  [£(n)]_1  •  Also, 
limn_>00[^(n)]~1  =  1 ,  converging  to  1  from  below  faster  than  (1  —  21~n). 

Thus,  as  n  grows  it  quickly  becomes  very  likely  that  n  randomly  selected  integers  have  a  gcd  of  1.  This 
fact,  together  with  Proposition  1,  make  estimation  of  t  via  our  algorithm  possible. 

The  parameter  estimation  techniques  given  above  lead  to  an  effective  method  for  periodic  pulse 
interval  analysis  (see  [CS13,  SC00,  SC98,  SC96,  Cas97,  CS96,  CS96]).  We  assume  time  is  highly  resolved 
and  ignore  any  time  quantization  error.  We  are  primarily  concerned  with  a  single  periodic  pulse  train 
with  (perhaps  very  many)  missing  observations  that  may  be  contaminated  with  outliers.  Our  data  model 
for  this  case,  in  terms  of  the  arrival  times  tj ,  is  given  as  above,  with  the  additional  assumption  that 
r)j  is  zero-mean  additive  white  Gaussian  noise.  Outliers  are  included  as  arbitrary  arrival  times.  The 
problem,  again,  is  to  recover  the  period  t  and  possibly  the  phase  <f.  With  Gaussian  noise  the  minimum 
variance  unbiased  estimates  for  this  linear  regression  problem  take  a  least-squares  form.  However,  this 
requires  knowledge  of  the  kj's.  We  therefore  have  developed  a  multi-step  procedure  that  proceeds  by  (i) 
estimating  r  directly,  (ii)  estimating  the  kj's,  and  (iii)  refining  the  estimate  of  r  using  the  estimated  kj's 
in  the  least-squares  solution. 

7.2  Multiple  Generators:  The  EQUIMEA  Algorithm 

We  close  by  discussing  our  work  on  deinterleaving.  Our  data  model  is  the  union  of  M  copies  of  our  previous 
datasets,  each  with  different  periods  or  “generators”  T  =  {r*},  kij's  and  phases.  Let  t  =  maxjjrj}.  Then 
our  data  is  S  =  U-=1  Wi  +  kijTi  +  ,  where  n,  is  the  number  of  elements  from  the  ith  generator, 

{kij}  is  a  linearly  increasing  sequence  of  natural  numbers  with  missing  observations,  <f>i  is  a  random 
variable  uniformly  distributed  in  [0,  Tj),  and  the  rjij's  are  zero-mean  iid  Gaussian  with  standard  deviation 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


New  Techniques  in  Time-Frequency  Analysis 


24 


3<7jj  <  t/2.  We  think  of  the  data  as  events  from  M  periodic  processes,  and  represent  it,  after  reindexing, 
as  S  =  {cq}^.  We  difference  as  in  the  MEA,  but  we  compute  all  of  the  differences.  We  repeat  this 
m  times,  and  saving  the  elements  from  each  iteration.  We  form  a  union  of  all  of  these  data  elements. 
The  relative  prinreness  of  data  generated  by  one  generator  will  “fill  in”  the  missing  elements  for  that 
generator,  whereas  the  data  from  two  different  generators  will  become  “Weyl  flat.”  The  number  of  times 
m  that  this  process  is  applied  will  be  determined  experimentally.  Assuming  only  minimal  knowledge  of 
the  range  of  {Tj},  namely  bounds  Tjj.  Ty  such  that  0  <Tjj  <  rt  <  Ty ,  we  phase  wrap  the  data  by  the 
mapping  4>p(az)  =  (y)  =  f  ~  f 
fractional  part,  and  so  4>p(cq)  E  [0, 1 


,  where  p  E  [Tl,Tjj\,  and  [-J  is  the  floor  function.  Thus  (■)  is  the 


Definition  9  A  sequence  of  real  random  variables  {xj}  C  [0, 1)  is  essentially  uniformly  distributed  in 
the  sense  of  Weyl  if  given  a,b,  0  <  a  <  b  <  1,  ^cardjl  <  j  <  n  :  Xj  E  [a,  6]}  — >  ( b  —  a)  as  n  — »  oo 
almost  surely. 


Weyl’s  Theorem  is  presented  in  [  ].  For  our  variation,  we  assume  that  for  each  i,  {kij}  is  a  linearly 

increasing  infinite  sequence  of  natural  numbers  with  missing  observations  such  that  kjj  — »  oo  as  j  — >  oo. 
We  must  make  this  assumption  because  the  result  is  only  approximately  true  for  a  finite  length  sequence. 


Theorem  15  For  almost  every  choice  of  p  (in  the  sense  of  Lebesgue  measure)  <&p(ai)  is  essentially 
uniformly  distributed  in  the  sense  of  Weyl. 

Moreover,  the  set  of  p’s  for  which  this  is  not  true  are  rational  multiples  of  {t*}.  Therefore,  except  for  those 
values,  is  essentially  uniformly  distributed  in  [0,1).  Moreover,  the  values  at  which  <&p(otij)  =  0 

almost  surely  are  p  E  {rj/n  :  n  E  N}.  These  values  of  p  cluster  at  zero,  but  spread  out  for  lower 
values  of  n.  We  phase  wrap  the  data  by  computing  modulus  of  the  spectrum,  i.e.,  compute  \Siter(r)\  = 
|  l  e^2ms^/T^\  .  The  values  of  |5'iter(r)|  will  have  peaks  at  the  periods  Tj  and  their  harmonics  (■ Tj)/k . 
The  “noise-like”  behavior  of  <hp(cq)  for  a.e.  p  leads  to  a  “flat”  range  for  S  for  p  0  {rj/n  :  n  E  N}.  In  turn, 
this  gives  the  following.  Let  io  denote  the  index  of  the  most  prolific  generator.  We  then  isolate  the  data 
generated  by  Tj0  by  convolution  with  a  pulse  train  of  width  Tj0,  and  subtract  it  out  of  Ao-  We  then  repeat 
the  process,  terminating  when  Ao  equals  the  empty  set.  We  refer  to  this  at  the  EQIMEA  algorithm. 

The  EQUIMEA  Algorithm  —  Multiple  Periods 


Initialize:  Sort  the  elements  of  S  in  descending  order.  Form  the  new  set  with  elements  (s;  —  sz+i)-  Set 
si  < —  (si  —  s;_|_i).  (Note,  this  eliminates  the  phase  ip.)  Set  iter  =  1,  i  =  1,  and  Error.  Go  to  1.) 

1. )  [Adjoin  0  after  first  iteration.]  Siter  * —  S  U  {0}. 

2. ) [Sort.]  Sort  the  elements  of  Siter  in  descending  order. 

3. ) [Compute  all  differences.]  Set  Siter  =  U(si  —  sk )  with  Sj  >  Sk- 

4. )  [Eliminate  zero(s).]  If  Sj  =  0,  then  Siter  * —  Suer  \  { s ^ - 

5. ) [Adjoin  previous  iteration.]  Form  Siter  * —  Siter  U  Siter- 1- 

6. )] [Compute  spectrum.]  Compute  \Siter{j)\  =  \  i  e^2ms^'>^T',  \ . 

7. )  [Threshold.]  Choose  the  largest  peak.  Label  it  as  Titer 

8. )  If  |  Tuer  —  Titer- 1|  <  Error.  Declare  fj  =  Titer ■  If  not,  iter  < —  (iter  +  1).  Go  to  1.). 

9. )  Given  Tj,  frequency  notch  | *S'?;^er (t) |  for  fi/m,  m  E  N.  Let  i  < —  i  +  1. 

10. )[Compute  spectrum.]  Compute  \Sit.er(T)\  =  \  J2f=i  \  . 

11. )  [Threshold.]  Choose  the  largest  peak.  Label  it  as  \.  Algorithm  terminates  when  there  are  no 
peaks.  Else,  go  to  9.). 
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Pulse  Train:  3  periods,  10%  of  data 


Figure  1:  Three  Periods  -  OriginalData 


We  now  demonstrate  the  algorithm.  The  original  data  had  three  underlying  periods  -1,(1  +  y/b)/2,  y/7, 
with  90%  of  the  information  randomly  removed,  and  10%  jitter  noise. 

Here  is  the  |Sl^er|  after  two  iterations. 


Figure  2:  EQUIMEA  Three  Periods  -  Iter2  -  Spectrum 

There  are  several  research  items  that  need  to  be  addressed.  First,  we  need  to  explore  applications 
of  both  algorithms  to  general  point  processes.  The  second  is  to  combine  this  approach  with  a  localized 
time-frequency  to  explore  signal  separation  algorithms.  We  plan  on  looking  at  waveform  signatures  in, 
e.g.,  EEG  data  (using  wavelets),  and  see  the  application  of  the  combination  of  time- frequency  analysis 
with  these  point  process  algorithms  to  the  deinterleaving  and  analysis. 
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